diff --git a/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h b/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h
index 8c8cf30dc0bd8..8907b200b8f3c 100644
--- a/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h
+++ b/roofit/histfactory/inc/RooStats/HistFactory/Measurement.h
@@ -650,9 +650,9 @@ extern Channel BadChannel;
 class Measurement : public TNamed {
 
 public:
-   Measurement();
-   ///  Measurement( const Measurement& other ); // Copy
-   Measurement(const char *Name, const char *Title = "");
+   Measurement() = default;
+   /// Constructor specifying name and title of measurement
+   Measurement(const char *Name, const char *Title = "") : TNamed(Name, Title) {}
 
    ///  set output prefix
    void SetOutputFilePrefix(const std::string &prefix) { fOutputFilePrefix = prefix; }
@@ -724,7 +724,8 @@ class Measurement : public TNamed {
    std::vector<RooStats::HistFactory::Channel> &GetChannels() { return fChannels; }
    const std::vector<RooStats::HistFactory::Channel> &GetChannels() const { return fChannels; }
    RooStats::HistFactory::Channel &GetChannel(std::string);
-   void AddChannel(RooStats::HistFactory::Channel chan);
+   /// Add a completely configured channel.
+   void AddChannel(RooStats::HistFactory::Channel chan) { fChannels.push_back(chan); }
 
    bool HasChannel(std::string);
    void writeToFile(TFile *file);
@@ -752,10 +753,10 @@ class Measurement : public TNamed {
    /// Configurables of this measurement
    std::string fOutputFilePrefix;
    std::vector<std::string> fPOI;
-   double fLumi;
-   double fLumiRelErr;
-   int fBinLow;
-   int fBinHigh;
+   double fLumi = 1.0;
+   double fLumiRelErr = 0.1;
+   int fBinLow = 0;
+   int fBinHigh = 1;
    bool fExportOnly = true;
    std::string fInterpolationScheme;
 
diff --git a/roofit/histfactory/src/Measurement.cxx b/roofit/histfactory/src/Measurement.cxx
index bc7b84e43e673..9955131003768 100644
--- a/roofit/histfactory/src/Measurement.cxx
+++ b/roofit/histfactory/src/Measurement.cxx
@@ -41,29 +41,6 @@ uncertainty or the functional form of constraints on nuisance parameters.
 
 namespace RooStats::HistFactory {
 
-/// Standard constructor
-Measurement::Measurement() : fLumi(1.0), fLumiRelErr(.10), fBinLow(0), fBinHigh(1) {}
-
-/*
-Measurement::Measurement(const Measurement& other) :
-  POI( other.POI ), Lumi( other.Lumi ), LumiRelErr( other.LumiRelErr ),
-  BinLow( other.BinLow ), BinHigh( other.BinHigh ), ExportOnly( other.ExportOnly ),
-  channels( other.channels ), OutputFilePrefix( other.outputFilePrefix ),
-  constantParams( other.constantParams ), { ; }
-*/
-
-/// Standard constructor specifying name and title of measurement
-Measurement::Measurement(const char *Name, const char *Title)
-   : TNamed(Name, Title), fLumi(1.0), fLumiRelErr(.10), fBinLow(0), fBinHigh(1)
-{
-}
-
-/// add a completely configured channel
-void Measurement::AddChannel(Channel chan)
-{
-   fChannels.push_back(chan);
-}
-
 /// Set a parameter in the model to be constant.
 /// the parameter does not have to exist yet, the information will be used when
 /// the model is actually created.
@@ -183,13 +160,6 @@ Channel &Measurement::GetChannel(std::string ChanName)
    // return BadChannel;
 }
 
-/*
-  void Measurement::Print( Option_t* option ) const {
-  Measurement::Print( std::cout );
-  return;
-  }
-*/
-
 /// Print information about measurement object in tree-like structure to given stream
 void Measurement::PrintTree(std::ostream &stream)
 {
@@ -270,14 +240,7 @@ void Measurement::PrintXML(std::string directory, std::string newOutputPrefix)
    // Add the time
    xml << "<!--" << std::endl;
    xml << "This xml file created automatically on: " << std::endl;
-   /*
-     time_t t = time(0);   // get time now
-     struct tm * now = localtime( &t );
-     xml << (now->tm_year + 1900) << '-'
-         << (now->tm_mon + 1) << '-'
-         << now->tm_mday
-         << std::endl;
-   */
+
    // LM: use TTimeStamp
    TTimeStamp t;
    UInt_t year = 0;
@@ -301,12 +264,6 @@ void Measurement::PrintXML(std::string directory, std::string newOutputPrefix)
    for (unsigned int i = 0; i < fFunctionObjects.size(); ++i) {
       PreprocessFunction func = fFunctionObjects.at(i);
       func.PrintXML(xml);
-      /*
-      xml << "<Function Name=\"" << func.GetName() << "\" "
-     << "Expression=\""     << func.GetExpression() << "\" "
-     << "Dependents=\""     << func.GetDependents() << "\" "
-     << "/>" << std::endl;
-      */
    }
 
    xml << std::endl;
@@ -418,15 +375,13 @@ void Measurement::writeToFile(TFile *file)
    // Collect all histograms from file:
    // HistCollector collector;
 
-   for (unsigned int chanItr = 0; chanItr < outMeas.fChannels.size(); ++chanItr) {
-
+   for (Channel &channel : fChannels) {
       // Go to the main directory
       // in the file
       file->cd();
       file->Flush();
 
       // Get the name of the channel:
-      Channel &channel = outMeas.fChannels.at(chanItr);
       std::string chanName = channel.GetName();
 
       if (!channel.CheckHistograms()) {
@@ -461,39 +416,21 @@ void Measurement::writeToFile(TFile *file)
 
       channel.fData.writeToFile(OutputFileName, GetDirPath(dataDir));
 
-      /*
-      // Write the data file to this directory
-      TH1* hData = channel.data.GetHisto();
-      hData->Write();
-
-      // Set the location of the data
-      // in the output measurement
-
-      channel.data.InputFile = OutputFileName;
-      channel.data.HistoName = hData->GetName();
-      channel.data.HistoPath = GetDirPath( dataDir );
-      */
-
       // Loop over samples:
-
-      for (unsigned int sampItr = 0; sampItr < channel.GetSamples().size(); ++sampItr) {
-
-         Sample &sample = channel.GetSamples().at(sampItr);
-         std::string sampName = sample.GetName();
-
-         cxcoutPHF << "Writing sample: " << sampName << std::endl;
+      for (Sample &sample : channel.GetSamples()) {
+         cxcoutPHF << "Writing sample: " << sample.GetName() << std::endl;
 
          file->cd();
          chanDir->cd();
-         TDirectory *sampleDir = chanDir->mkdir(sampName.c_str());
+         TDirectory *sampleDir = chanDir->mkdir(sample.GetName().c_str());
          if (sampleDir == nullptr) {
-            cxcoutEHF << "Error: Directory " << sampName << " not created properly" << std::endl;
+            cxcoutEHF << "Error: Directory " << sample.GetName() << " not created properly" << std::endl;
             throw hf_exc();
          }
          std::string sampleDirPath = GetDirPath(sampleDir);
 
          if (!sampleDir) {
-            cxcoutEHF << "Error making directory: " << sampName << " in directory: " << chanName << std::endl;
+            cxcoutEHF << "Error making directory: " << sample.GetName() << " in directory: " << chanName << std::endl;
             throw hf_exc();
          }
 
@@ -501,59 +438,6 @@ void Measurement::writeToFile(TFile *file)
          sampleDir->cd();
 
          sample.writeToFile(OutputFileName, sampleDirPath);
-         /*
-         TH1* hSample = sample.GetHisto();
-         if( ! hSample ) {
-      std::cout << "Error getting histogram for sample: "
-           << sampName << std::endl;
-      throw -1;
-         }
-         sampleDir->cd();
-         hSample->Write();
-
-         sample.InputFile = OutputFileName;
-         sample.HistoName = hSample->GetName();
-         sample.HistoPath = sampleDirPath;
-         */
-
-         // Write the histograms associated with
-         // systematics
-
-         /*  THIS IS WHAT I"M COMMENTING
-         sample.GetStatError().writeToFile( OutputFileName, sampleDirPath );
-
-         // Must write all systematics that contain internal histograms
-         // (This is not all systematics)
-
-         for( unsigned int i = 0; i < sample.GetHistoSysList().size(); ++i ) {
-      sample.GetHistoSysList().at(i).writeToFile( OutputFileName, sampleDirPath );
-         }
-         for( unsigned int i = 0; i < sample.GetHistoFactorList().size(); ++i ) {
-      sample.GetHistoFactorList().at(i).writeToFile( OutputFileName, sampleDirPath );
-         }
-         for( unsigned int i = 0; i < sample.GetShapeSysList().size(); ++i ) {
-      sample.GetShapeSysList().at(i).writeToFile( OutputFileName, sampleDirPath );
-         }
-         END COMMENT  */
-         /*
-         sample.statError.writeToFile( OutputFileName, sampleDirPath );
-
-         // Now, get the Stat config histograms
-         if( sample.statError.HistoName != "" ) {
-      TH1* hStatError = sample.statError.GetErrorHist();
-      if( ! hStatError ) {
-        std::cout << "Error getting stat error histogram for sample: "
-             << sampName << std::endl;
-        throw -1;
-      }
-      hStatError->Write();
-
-      sample.statError.InputFile = OutputFileName;
-      sample.statError.HistoName = hStatError->GetName();
-      sample.statError.HistoPath = sampleDirPath;
-
-         }
-         */
       }
    }
 
@@ -571,7 +455,6 @@ void Measurement::writeToFile(TFile *file)
 /// stripped of unnecessary prefixes
 std::string Measurement::GetDirPath(TDirectory *dir)
 {
-
    std::string path = dir->GetPath();
 
    if (path.find(':') != std::string::npos) {
@@ -579,18 +462,7 @@ std::string Measurement::GetDirPath(TDirectory *dir)
       path.replace(0, index + 1, "");
    }
 
-   path = path + "/";
-
-   return path;
-
-   /*
-       if( path.find(":") != std::string::npos ) {
-    size_t index = path.find(":");
-    SampleName.replace( 0, index, "" );
-       }
-
-       // Remove the file:
-       */
+   return path + "/";
 }
 
 /// The most common way to add histograms to channels is to have them
@@ -603,11 +475,7 @@ std::string Measurement::GetDirPath(TDirectory *dir)
 /// save all necessary histograms.
 void Measurement::CollectHistograms()
 {
-
-   for (unsigned int chanItr = 0; chanItr < fChannels.size(); ++chanItr) {
-
-      Channel &chan = fChannels.at(chanItr);
-
+   for (Channel &chan : fChannels) {
       chan.CollectHistograms();
    }
 }
@@ -727,8 +595,6 @@ void Sample::writeToFile(std::string OutputFileName, std::string DirName)
    for (unsigned int i = 0; i < GetShapeFactorList().size(); ++i) {
       GetShapeFactorList().at(i).writeToFile(OutputFileName, DirName);
    }
-
-   return;
 }
 
 void Sample::SetValue(double val)
@@ -1200,10 +1066,7 @@ void Channel::CollectHistograms()
    }
 
    // Get the histograms for the samples:
-   for (unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr) {
-
-      Sample &sample = fSamples.at(sampItr);
-
+   for (Sample &sample : fSamples) {
       // Get the nominal histogram:
       cxcoutDHF << "Collecting Nominal Histogram" << std::endl;
       TH1 *Nominal = GetHistogram(sample.GetInputFile(), sample.GetHistoPath(), sample.GetHistoName(), fileHandles);
@@ -1218,10 +1081,7 @@ void Channel::CollectHistograms()
       }
 
       // Get the HistoSys Variations:
-      for (unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr) {
-
-         HistoSys &histoSys = sample.GetHistoSysList().at(histoSysItr);
-
+      for (HistoSys &histoSys : sample.GetHistoSysList()) {
          histoSys.SetHistoLow(GetHistogram(histoSys.GetInputFileLow(), histoSys.GetHistoPathLow(),
                                            histoSys.GetHistoNameLow(), fileHandles));
 
@@ -1230,10 +1090,7 @@ void Channel::CollectHistograms()
       } // End Loop over HistoSys
 
       // Get the HistoFactor Variations:
-      for (unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr) {
-
-         HistoFactor &histoFactor = sample.GetHistoFactorList().at(histoFactorItr);
-
+      for (HistoFactor &histoFactor : sample.GetHistoFactorList()) {
          histoFactor.SetHistoLow(GetHistogram(histoFactor.GetInputFileLow(), histoFactor.GetHistoPathLow(),
                                               histoFactor.GetHistoNameLow(), fileHandles));
 
@@ -1242,19 +1099,13 @@ void Channel::CollectHistograms()
       } // End Loop over HistoFactor
 
       // Get the ShapeSys Variations:
-      for (unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr) {
-
-         ShapeSys &shapeSys = sample.GetShapeSysList().at(shapeSysItr);
-
+      for (ShapeSys &shapeSys : sample.GetShapeSysList()) {
          shapeSys.SetErrorHist(
             GetHistogram(shapeSys.GetInputFile(), shapeSys.GetHistoPath(), shapeSys.GetHistoName(), fileHandles));
       } // End Loop over ShapeSys
 
       // Get any initial shape for a ShapeFactor
-      for (unsigned int shapeFactorItr = 0; shapeFactorItr < sample.GetShapeFactorList().size(); ++shapeFactorItr) {
-
-         ShapeFactor &shapeFactor = sample.GetShapeFactorList().at(shapeFactorItr);
-
+      for (ShapeFactor &shapeFactor : sample.GetShapeFactorList()) {
          // Check if we need an InitialShape
          if (shapeFactor.HasInitialShape()) {
             TH1 *hist = GetHistogram(shapeFactor.GetInputFile(), shapeFactor.GetHistoPath(), shapeFactor.GetHistoName(),
@@ -1279,10 +1130,7 @@ bool Channel::CheckHistograms() const
    }
 
    // Get the histograms for the samples:
-   for (unsigned int sampItr = 0; sampItr < fSamples.size(); ++sampItr) {
-
-      const Sample &sample = fSamples.at(sampItr);
-
+   for (Sample const &sample : fSamples) {
       // Get the nominal histogram:
       if (sample.GetHisto() == nullptr) {
          cxcoutEHF << "Error: Nominal Histogram for sample " << sample.GetName() << " is nullptr." << std::endl;
@@ -1322,9 +1170,7 @@ bool Channel::CheckHistograms() const
       }
 
       // Get the HistoSys Variations:
-      for (unsigned int histoSysItr = 0; histoSysItr < sample.GetHistoSysList().size(); ++histoSysItr) {
-
-         const HistoSys &histoSys = sample.GetHistoSysList().at(histoSysItr);
+      for (const HistoSys &histoSys : sample.GetHistoSysList()) {
 
          if (histoSys.GetHistoLow() == nullptr) {
             cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoSys.GetName() << " in sample "
@@ -1340,9 +1186,7 @@ bool Channel::CheckHistograms() const
       } // End Loop over HistoSys
 
       // Get the HistoFactor Variations:
-      for (unsigned int histoFactorItr = 0; histoFactorItr < sample.GetHistoFactorList().size(); ++histoFactorItr) {
-
-         const HistoFactor &histoFactor = sample.GetHistoFactorList().at(histoFactorItr);
+      for (const HistoFactor &histoFactor : sample.GetHistoFactorList()) {
 
          if (histoFactor.GetHistoLow() == nullptr) {
             cxcoutEHF << "Error: HistoSyst Low for Systematic " << histoFactor.GetName() << " in sample "
@@ -1358,16 +1202,12 @@ bool Channel::CheckHistograms() const
       } // End Loop over HistoFactor
 
       // Get the ShapeSys Variations:
-      for (unsigned int shapeSysItr = 0; shapeSysItr < sample.GetShapeSysList().size(); ++shapeSysItr) {
-
-         const ShapeSys &shapeSys = sample.GetShapeSysList().at(shapeSysItr);
-
+      for (const ShapeSys &shapeSys : sample.GetShapeSysList()) {
          if (shapeSys.GetErrorHist() == nullptr) {
             cxcoutEHF << "Error: HistoSyst High for Systematic " << shapeSys.GetName() << " in sample "
                       << sample.GetName() << " is nullptr." << std::endl;
             return false;
          }
-
       } // End Loop over ShapeSys
 
    } // End Loop over Samples
@@ -1614,10 +1454,7 @@ void Asimov::ConfigureWorkspace(RooWorkspace *wspace)
    // Then, we set any variables to constant
    //
 
-   for (std::map<std::string, bool>::iterator itr = fParamsToFix.begin(); itr != fParamsToFix.end(); ++itr) {
-
-      std::string param = itr->first;
-      bool isConstant = itr->second;
+   for (auto &[param, isConstant] : fParamsToFix) {
 
       // Try to get the variable in the workspace
       RooRealVar *var = wspace->var(param);
@@ -1630,8 +1467,6 @@ void Asimov::ConfigureWorkspace(RooWorkspace *wspace)
       std::cout << "Configuring Asimov Dataset: Setting " << param << " to constant " << std::endl;
       var->setConstant(isConstant);
    }
-
-   return;
 }
 
 /** \class RooStats::HistFactory::HistRef
@@ -1737,8 +1572,6 @@ Constraint::Type Constraint::GetType(const std::string &Name)
    }
 }
 
-// Norm Factor
-
 void NormFactor::Print(std::ostream &stream) const
 {
    stream << "\t \t Name: " << fName << "\t Val: " << fVal << "\t Low: " << fLow << "\t High: " << fHigh << std::endl;
@@ -1753,7 +1586,6 @@ void NormFactor::PrintXML(std::ostream &xml) const
        << "  /> " << std::endl;
 }
 
-// Overall Sys
 void OverallSys::Print(std::ostream &stream) const
 {
    stream << "\t \t Name: " << fName << "\t Low: " << fLow << "\t High: " << fHigh << std::endl;
@@ -1814,8 +1646,6 @@ void HistogramUncertaintyBase::writeToFile(const std::string &FileName, const st
    fInputFileHigh = FileName;
    fHistoPathHigh = DirName;
    fHistoNameHigh = histHigh->GetName();
-
-   return;
 }
 
 void HistoSys::PrintXML(std::ostream &xml) const
@@ -1831,8 +1661,6 @@ void HistoSys::PrintXML(std::ostream &xml) const
        << "  /> " << std::endl;
 }
 
-// Shape Sys
-
 void ShapeSys::SetErrorHist(TH1 *hError)
 {
    fhHigh.reset(hError);
@@ -1856,7 +1684,6 @@ void ShapeSys::PrintXML(std::ostream &xml) const
 
 void ShapeSys::writeToFile(const std::string &FileName, const std::string &DirName)
 {
-
    auto histError = GetErrorHist();
    if (histError == nullptr) {
       std::cout << "Error: Cannot write " << GetName() << " to file: " << FileName << " ErrorHist is nullptr"
@@ -1867,12 +1694,8 @@ void ShapeSys::writeToFile(const std::string &FileName, const std::string &DirNa
    fInputFileHigh = FileName;
    fHistoPathHigh = DirName;
    fHistoNameHigh = histError->GetName();
-
-   return;
 }
 
-// HistoFactor
-
 void HistoFactor::PrintXML(std::ostream &xml) const
 {
    xml << "      <HistoFactor Name=\"" << GetName() << "\" "
@@ -1887,8 +1710,6 @@ void HistoFactor::PrintXML(std::ostream &xml) const
        << "  /> " << std::endl;
 }
 
-// Shape Factor
-
 void ShapeFactor::SetInitialShape(TH1 *shape)
 {
    fhHigh.reset(shape);
@@ -1925,8 +1746,6 @@ void ShapeFactor::writeToFile(const std::string &FileName, const std::string &Di
       fHistoPathHigh = DirName;
       fHistoNameHigh = histInitialShape->GetName();
    }
-
-   return;
 }
 
 void ShapeFactor::PrintXML(std::ostream &xml) const
@@ -1945,7 +1764,6 @@ void StatError::SetErrorHist(TH1 *Error)
    fhHigh.reset(Error);
 }
 
-// Stat Error Config
 void StatErrorConfig::Print(std::ostream &stream) const
 {
    stream << "\t \t RelErrorThreshold: " << fRelErrorThreshold
@@ -1960,7 +1778,6 @@ void StatErrorConfig::PrintXML(std::ostream &xml) const
        << std::endl;
 }
 
-// Stat Error
 void StatError::Print(std::ostream &stream) const
 {
    stream << "\t \t Activate: " << fActivate << "\t InputFile: " << fInputFileHigh << "\t HistoName: " << fHistoNameHigh
@@ -1997,8 +1814,6 @@ void StatError::writeToFile(const std::string &OutputFileName, const std::string
       fHistoNameHigh = statErrorHistName;
       fHistoPathHigh = DirName;
    }
-
-   return;
 }
 
 HistogramUncertaintyBase::HistogramUncertaintyBase(const HistogramUncertaintyBase &oth)
diff --git a/roofit/roofitcore/inc/RooFitResult.h b/roofit/roofitcore/inc/RooFitResult.h
index 2d2181ee4adec..0eabae7df175e 100644
--- a/roofit/roofitcore/inc/RooFitResult.h
+++ b/roofit/roofitcore/inc/RooFitResult.h
@@ -33,9 +33,7 @@
 class RooArgSet ;
 class RooAbsPdf ;
 class RooPlot;
-class TObject ;
-class TH2 ;
-typedef RooArgSet* pRooArgSet ;
+class TH2;
 
 class RooFitResult : public TNamed, public RooPrintable, public RooDirItem {
 public:
@@ -185,11 +183,12 @@ class RooFitResult : public TNamed, public RooPrintable, public RooDirItem {
   RooArgList* _initPars = nullptr;     ///< List of floating parameters with initial values
   RooArgList* _finalPars = nullptr;    ///< List of floating parameters with final values
 
-  mutable RooArgList* _globalCorr = nullptr;   ///<! List of global correlation coefficients
+  mutable std::unique_ptr<RooArgList> _globalCorr; ///<! List of global correlation coefficients
   mutable TList       _corrMatrix ;   ///<! Correlation matrix (list of RooArgLists)
 
-  mutable RooArgList *_randomPars = nullptr; ///<! List of floating parameters with most recent random perturbation applied
-  mutable TMatrixF* _Lt = nullptr;           ///<! triangular matrix used for generate random perturbations
+  mutable std::unique_ptr<RooArgList>
+     _randomPars;                        ///<! List of floating parameters with most recent random perturbation applied
+  mutable std::unique_ptr<TMatrixF> _Lt; ///<! triangular matrix used for generate random perturbations
 
   TMatrixDSym* _CM = nullptr;  ///< Correlation matrix
   TMatrixDSym* _VM = nullptr;  ///< Covariance matrix
diff --git a/roofit/roofitcore/src/RooFitResult.cxx b/roofit/roofitcore/src/RooFitResult.cxx
index b86417a372c6c..a1900eb46f6aa 100644
--- a/roofit/roofitcore/src/RooFitResult.cxx
+++ b/roofit/roofitcore/src/RooFitResult.cxx
@@ -92,10 +92,11 @@ RooFitResult::RooFitResult(const RooFitResult &other)
 
   other._finalPars->snapshot(*_finalPars);
   if (other._randomPars) {
-    _randomPars = new RooArgList;
-    other._randomPars->snapshot(*_randomPars);
+     _randomPars = std::make_unique<RooArgList>();
+     other._randomPars->snapshot(*_randomPars);
   }
-  if (other._Lt) _Lt = new TMatrix(*other._Lt);
+  if (other._Lt)
+     _Lt = std::make_unique<TMatrix>(*other._Lt);
   if (other._VM) _VM = new TMatrixDSym(*other._VM) ;
   if (other._CM) _CM = new TMatrixDSym(*other._CM) ;
   if (other._GC) _GC = new TVectorD(*other._GC) ;
@@ -113,10 +114,8 @@ RooFitResult::~RooFitResult()
 {
   if (_constPars) delete _constPars ;
   if (_initPars)  delete _initPars ;
-  if (_finalPars) delete _finalPars ;
-  if (_globalCorr) delete _globalCorr;
-  if (_randomPars) delete _randomPars;
-  if (_Lt) delete _Lt;
+  if (_finalPars)
+     delete _finalPars;
   if (_CM) delete _CM ;
   if (_VM) delete _VM ;
   if (_GC) delete _GC ;
@@ -334,37 +333,36 @@ RooPlot *RooFitResult::plotOn(RooPlot *frame, const char *parName1, const char *
 const RooArgList& RooFitResult::randomizePars() const
 {
   Int_t nPar= _finalPars->size();
-  if(nullptr == _randomPars) { // first-time initialization
-    assert(nullptr != _finalPars);
-    // create the list of random values to fill
-    _randomPars = new RooArgList;
-    _finalPars->snapshot(*_randomPars);
-    // calculate the elements of the upper-triangular matrix L that gives Lt*L = C
-    // where Lt is the transpose of L (the "square-root method")
-    TMatrix L(nPar,nPar);
-    for(Int_t iPar= 0; iPar < nPar; iPar++) {
-      // calculate the diagonal term first
-      L(iPar,iPar)= covariance(iPar,iPar);
-      for(Int_t k= 0; k < iPar; k++) {
-   double tmp= L(k,iPar);
-   L(iPar,iPar)-= tmp*tmp;
-      }
-      L(iPar,iPar)= sqrt(L(iPar,iPar));
-      // then the off-diagonal terms
-      for(Int_t jPar= iPar+1; jPar < nPar; jPar++) {
-   L(iPar,jPar)= covariance(iPar,jPar);
-   for(Int_t k= 0; k < iPar; k++) {
-     L(iPar,jPar)-= L(k,iPar)*L(k,jPar);
-   }
-   L(iPar,jPar)/= L(iPar,iPar);
-      }
-    }
-    // remember Lt
-    _Lt= new TMatrix(TMatrix::kTransposed,L);
-  }
-  else {
-    // reset to the final fit values
-    _randomPars->assign(*_finalPars);
+  if (!_randomPars) { // first-time initialization
+     assert(nullptr != _finalPars);
+     // create the list of random values to fill
+     _randomPars = std::make_unique<RooArgList>();
+     _finalPars->snapshot(*_randomPars);
+     // calculate the elements of the upper-triangular matrix L that gives Lt*L = C
+     // where Lt is the transpose of L (the "square-root method")
+     TMatrix L(nPar, nPar);
+     for (Int_t iPar = 0; iPar < nPar; iPar++) {
+        // calculate the diagonal term first
+        L(iPar, iPar) = covariance(iPar, iPar);
+        for (Int_t k = 0; k < iPar; k++) {
+           double tmp = L(k, iPar);
+           L(iPar, iPar) -= tmp * tmp;
+        }
+        L(iPar, iPar) = sqrt(L(iPar, iPar));
+        // then the off-diagonal terms
+        for (Int_t jPar = iPar + 1; jPar < nPar; jPar++) {
+           L(iPar, jPar) = covariance(iPar, jPar);
+           for (Int_t k = 0; k < iPar; k++) {
+              L(iPar, jPar) -= L(k, iPar) * L(k, jPar);
+           }
+           L(iPar, jPar) /= L(iPar, iPar);
+        }
+     }
+     // remember Lt
+     _Lt = std::make_unique<TMatrix>(TMatrix::kTransposed, L);
+  } else {
+     // reset to the final fit values
+     _randomPars->assign(*_finalPars);
   }
 
   // create a vector of unit Gaussian variables
@@ -455,7 +453,7 @@ const RooArgList* RooFitResult::globalCorr()
     fillLegacyCorrMatrix() ;
   }
 
-  return _globalCorr ;
+  return _globalCorr.get();
 }
 
 
@@ -622,39 +620,28 @@ void RooFitResult::fillLegacyCorrMatrix() const
   if (!_CM) return ;
 
   // Delete eventual previous correlation data holders
-  if (_globalCorr) delete _globalCorr ;
   _corrMatrix.Delete();
 
   // Build holding arrays for correlation coefficients
-  _globalCorr = new RooArgList("globalCorrelations") ;
+  _globalCorr = std::make_unique<RooArgList>("globalCorrelations");
 
   for(RooAbsArg* arg : *_initPars) {
     // Create global correlation value holder
-    TString gcName("GC[") ;
-    gcName.Append(arg->GetName()) ;
-    gcName.Append("]") ;
-    TString gcTitle(arg->GetTitle()) ;
-    gcTitle.Append(" Global Correlation") ;
-    _globalCorr->addOwned(std::make_unique<RooRealVar>(gcName.Data(),gcTitle.Data(),0.));
+    std::string argName = arg->GetName();
+    std::string argTitle = arg->GetTitle();
+    std::string gcName = "GC[" + argName + "]";
+    std::string gcTitle = argTitle + " Global Correlation";
+    _globalCorr->addOwned(std::make_unique<RooRealVar>(gcName.c_str(), gcTitle.c_str(), 0.));
 
     // Create array with correlation holders for this parameter
-    TString name("C[") ;
-    name.Append(arg->GetName()) ;
-    name.Append(",*]") ;
-    RooArgList* corrMatrixRow = new RooArgList(name.Data()) ;
+    RooArgList *corrMatrixRow = new RooArgList(("C[" + argName + ",*]").c_str());
     _corrMatrix.Add(corrMatrixRow) ;
     for(RooAbsArg* arg2 : *_initPars) {
 
-      TString cName("C[") ;
-      cName.Append(arg->GetName()) ;
-      cName.Append(",") ;
-      cName.Append(arg2->GetName()) ;
-      cName.Append("]") ;
-      TString cTitle("Correlation between ") ;
-      cTitle.Append(arg->GetName()) ;
-      cTitle.Append(" and ") ;
-      cTitle.Append(arg2->GetName()) ;
-      corrMatrixRow->addOwned(std::make_unique<RooRealVar>(cName.Data(),cTitle.Data(),0.));
+       std::string arg2Name = arg2->GetName();
+       std::string cName = "C[" + argName + "," + arg2Name + "]";
+       std::string cTitle = "Correlation between " + argName + " and " + arg2Name;
+       corrMatrixRow->addOwned(std::make_unique<RooRealVar>(cName.c_str(), cTitle.c_str(), 0.));
     }
   }
 
@@ -1315,7 +1302,9 @@ void RooFitResult::Streamer(TBuffer &R__b)
       R__b >> _constPars;
       R__b >> _initPars;
       R__b >> _finalPars;
-      R__b >> _globalCorr;
+      RooArgList *globalCorr;
+      R__b >> globalCorr;
+      _globalCorr.reset(globalCorr);
       _corrMatrix.Streamer(R__b);
       R__b.CheckByteCount(R__s, R__c, RooFitResult::IsA());
 
diff --git a/roofit/roofitcore/src/RooPlot.cxx b/roofit/roofitcore/src/RooPlot.cxx
index 89cf8818476eb..a3a833b2e7cb0 100644
--- a/roofit/roofitcore/src/RooPlot.cxx
+++ b/roofit/roofitcore/src/RooPlot.cxx
@@ -898,22 +898,28 @@ bool RooPlot::drawAfter(const char *after, const char *target)
 /// methods to change the drawing style attributes of a contained
 /// object directly.
 
-TObject *RooPlot::findObject(const char *name, const TClass* tClass) const
+TObject *RooPlot::findObject(const char *name, const TClass *tClass) const
 {
-  TObject *ret = nullptr;
+   TObject *ret = nullptr;
 
-  for(auto const& item : _items) {
-    TObject &obj = *item.first;
-    if ((!name || name[0] == '\0' || !TString(name).CompareTo(obj.GetName()))
-        && (!tClass || (obj.IsA()==tClass))) {
-      ret = &obj ;
-    }
-  }
+   for (auto const &item : _items) {
+      TObject &obj = *item.first;
+      if ((!name || name[0] == '\0' || !TString(name).CompareTo(obj.GetName())) && (!tClass || (obj.IsA() == tClass))) {
+         ret = &obj;
+      }
+   }
 
-  if (ret == nullptr) {
-    coutE(InputArguments) << "RooPlot::findObject(" << GetName() << ") cannot find object " << (name?name:"<last>") << std::endl ;
-  }
-  return ret ;
+   if (ret == nullptr) {
+      std::stringstream error;
+      error << "RooPlot::findObject(" << GetName() << ") cannot find object " << (name ? name : "<last>") << "\n"
+            << "Available objects are:\n";
+      for (auto const &item : _items) {
+         TObject &obj = *item.first;
+         error << "  - " << obj.IsA()->GetName() << " \"" << obj.GetName() << "\"\n";
+      }
+      coutE(InputArguments) << error.str();
+   }
+   return ret;
 }
 
 
diff --git a/roofit/roofitcore/src/RooSimultaneous.cxx b/roofit/roofitcore/src/RooSimultaneous.cxx
index 0b241f961c493..7169024564fde 100644
--- a/roofit/roofitcore/src/RooSimultaneous.cxx
+++ b/roofit/roofitcore/src/RooSimultaneous.cxx
@@ -777,7 +777,8 @@ RooPlot* RooSimultaneous::plotOn(RooPlot *frame, RooLinkedList& cmdList) const
        projDataTmp(projData->reduce(RooFit::SelectVars(projDataVars), RooFit::Cut(cutString.c_str())));
 
     // Override normalization and projection dataset
-    RooCmdArg tmp1 = RooFit::Normalization(scaleFactor*wTable->getFrac(idxCatClone->getCurrentLabel()),stype) ;
+    RooCmdArg tmp1 =
+       RooFit::Normalization(scaleFactor * wTable->get(idxCatClone->getCurrentLabel()), RooAbsReal::NumEvent);
     RooCmdArg tmp2 = RooFit::ProjWData(*projDataSet,*projDataTmp) ;
 
     // WVE -- do not adjust normalization for asymmetry plots
diff --git a/roofit/roofitcore/test/gtest_wrapper.h b/roofit/roofitcore/test/gtest_wrapper.h
index de31e562ad9f5..f342c55575ab8 100644
--- a/roofit/roofitcore/test/gtest_wrapper.h
+++ b/roofit/roofitcore/test/gtest_wrapper.h
@@ -32,6 +32,28 @@
 
 #define ROOFIT_EVAL_BACKENDS ROOFIT_EVAL_BACKEND_LEGACY ROOFIT_EVAL_BACKEND_CUDA RooFit::EvalBackend::Cpu()
 
-#define ROOFIT_EVAL_BACKENDS_WITH_CODEGEN ROOFIT_EVAL_BACKENDS, ROOFIT_EVAL_BACKEND_CODEGEN RooFit::EvalBackend::CodegenNoGrad()
+#define ROOFIT_EVAL_BACKENDS_WITH_CODEGEN \
+   ROOFIT_EVAL_BACKENDS, ROOFIT_EVAL_BACKEND_CODEGEN RooFit::EvalBackend::CodegenNoGrad()
+
+#include <gtest/gtest.h>
+#include <gmock/gmock.h>
+
+#include <cmath>
+
+MATCHER_P2(RelativeNear, expected, rel_tol,
+           "is within relative tolerance around ref=" + ::testing::PrintToString(expected) +
+              " (tol=" + ::testing::PrintToString(rel_tol) + ")")
+{
+   const double diff = std::fabs(arg - expected);
+   const double scale = std::fabs(expected);
+   // One could also also consider the target value as scale, or an adaptive scale:
+   // const double scale = std::max(std::fabs(arg), std::fabs(expected));
+
+   if (diff <= rel_tol * scale)
+      return true;
+
+   *result_listener << "error relative to ref = " << diff / scale << ", absolute diff = " << diff;
+   return false;
+}
 
 #endif // RooFit_gtest_wrapper_h
diff --git a/roofit/roofitcore/test/stressRooFit_tests.h b/roofit/roofitcore/test/stressRooFit_tests.h
index 652414c05acd3..6cb169eedf035 100644
--- a/roofit/roofitcore/test/stressRooFit_tests.h
+++ b/roofit/roofitcore/test/stressRooFit_tests.h
@@ -2972,22 +2972,21 @@ class TestBasic501 : public RooUnitTest {
       RooPlot *frame1 = x.frame(Bins(30), Title("Physics sample"));
 
       // Plot all data tagged as physics sample
-      combData.plotOn(frame1, Cut("sample==sample::physics"));
+      std::unique_ptr<RooAbsData> slicedData1{combData.reduce(Cut("sample==sample::physics"))};
+      slicedData1->SetName("combData_Cut[sample==sample::physics]"); // to be consistent with reference file
+      slicedData1->plotOn(frame1);
 
       // Plot "physics" slice of simultaneous pdf.
-      // NBL You _must_ project the sample index category with data using ProjWData
-      // as a RooSimultaneous makes no prediction on the shape in the index category
-      // and can thus not be integrated
-      simPdf.plotOn(frame1, Slice(sample, "physics"), ProjWData(sample, combData));
-      simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData),
-                    LineStyle(kDashed));
-
-      // The same plot for the control sample slice
+      simPdf.getPdf("physics")->plotOn(frame1);
+      simPdf.getPdf("physics")->plotOn(frame1, Components("px"), LineStyle(kDashed));
+
+      // The same plot for the control sample slice.
       RooPlot *frame2 = x.frame(Bins(30), Title("Control sample"));
-      combData.plotOn(frame2, Cut("sample==sample::control"));
-      simPdf.plotOn(frame2, Slice(sample, "control"), ProjWData(sample, combData));
-      simPdf.plotOn(frame2, Slice(sample, "control"), Components("px_ctl"), ProjWData(sample, combData),
-                    LineStyle(kDashed));
+      std::unique_ptr<RooAbsData> slicedData2{combData.reduce(Cut("sample==sample::control"))};
+      slicedData2->SetName("combData_Cut[sample==sample::control]"); // to be consistent with reference file
+      slicedData2->plotOn(frame2);
+      simPdf.getPdf("control")->plotOn(frame2);
+      simPdf.getPdf("control")->plotOn(frame2, Components("px_ctl"), LineStyle(kDashed));
 
       regPlot(frame1, "rf501_plot1");
       regPlot(frame2, "rf501_plot2");
diff --git a/roofit/roofitcore/test/testRooSimultaneous.cxx b/roofit/roofitcore/test/testRooSimultaneous.cxx
index 14ac5e25ac074..548d695dc3ff3 100644
--- a/roofit/roofitcore/test/testRooSimultaneous.cxx
+++ b/roofit/roofitcore/test/testRooSimultaneous.cxx
@@ -5,10 +5,12 @@
 #include <RooAddPdf.h>
 #include <RooAddition.h>
 #include <RooCategory.h>
+#include <RooChebychev.h>
 #include <RooConstVar.h>
 #include <RooDataSet.h>
 #include <RooExponential.h>
 #include <RooFitResult.h>
+#include <RooGaussian.h>
 #include <RooGenericPdf.h>
 #include <RooHelpers.h>
 #include <RooMinimizer.h>
@@ -417,6 +419,8 @@ TEST(RooSimultaneous, ConditionalProdPdf)
 // channels can be extended. Also check if the likelihood can be created.
 TEST(RooSimultaneous, PartiallyExtendedPdfs)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
+
    RooWorkspace ws;
    ws.factory("Gaussian::pdfA(x_a[-10, 10], mu_a[0, -10, 10], sigma_a[2.0, 0.1, 10.0])");
    ws.factory("Gaussian::pdfB(x_b[-10, 10], mu_b[0, -10, 10], sigma_b[2.0, 0.1, 10.0])");
@@ -427,7 +431,6 @@ TEST(RooSimultaneous, PartiallyExtendedPdfs)
    RooArgSet observables{*ws.var("x_a"), *ws.var("x_b"), *ws.cat("cat")};
 
    auto &simPdf = *ws.pdf("simPdf");
-   std::cout << simPdf.getVal() << std::endl;
 
    // A completely extended pdf, just to easily create a toy dataset
    ws.factory("ExtendPdf::pdfAext(pdfA, n_b[1000., 100., 10000.])");
@@ -518,6 +521,7 @@ TEST_P(TestStatisticTest, RooSimultaneousSingleChannelCrossCheckWithCondVar)
 /// RooSimultaneous with the new CPU backend.
 TEST(RooSimultaneous, RangedExtendedRooAddPdf)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
 
    const double nBkgA_nom = 9000;
    const double nBkgB_nom = 10000;
@@ -549,11 +553,12 @@ TEST(RooSimultaneous, RangedExtendedRooAddPdf)
    simPdf.getParameters(combData->get(), params);
    params.snapshot(paramsSnap);
 
-   Res fitResult{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(RooFit::EvalBackend::Cpu()))};
+   Res fitResult{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Cpu()), PrintLevel(-1))};
 
    params.assign(paramsSnap);
 
-   Res fitResultRef{simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(RooFit::EvalBackend::Legacy()))};
+   Res fitResultRef{
+      simPdf.fitTo(*combData, Save(), Range("fitRange"), EvalBackend(EvalBackend::Legacy()), PrintLevel(-1))};
 
    EXPECT_TRUE(fitResult->isIdentical(*fitResultRef));
 }
@@ -563,6 +568,8 @@ TEST(RooSimultaneous, RangedExtendedRooAddPdf)
 /// with a projection dataset.
 TEST(RooSimultaneous, PlotProjWData)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLvl(RooFit::WARNING);
+
    RooRealVar x("x", "x", -8, 8);
    x.setBins(1);
 
@@ -590,11 +597,147 @@ TEST(RooSimultaneous, PlotProjWData)
    EXPECT_DOUBLE_EQ(frame->getCurve()->interpolate(0.), combData.sumEntries());
 }
 
+/// Second part of GitHub issue #20383.
+/// Check that the the simultaneous pdf is normalized correctly to the data
+/// when plotting with a projection dataset, in the extended and non-extended
+/// case, based on the reproducer provided by the user who opened the issue.
+TEST(RooSimultaneous, PlotProjWDataExtended)
+{
+   using namespace RooFit;
+
+   RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
+
+   RooRealVar xvar1("x1", "", 0.0, 10.0);
+   RooRealVar xvar2("x2", "", 0.0, 10.0);
+
+   RooRealVar mean1("mean1", "", 3.0, 1.0, 4.0);
+   RooRealVar mean2("mean2", "", 5.0, 4.0, 6.0);
+   RooRealVar sigma("sigma", "", 1.0, 0.01, 2.0);
+
+   RooGaussian gauss1("gauss1", "", xvar1, mean1, sigma);
+   RooGaussian gauss2("gauss2", "", xvar2, mean2, sigma);
+
+   RooRealVar a0("a0", "a0", -0.1, -1.0, 1.0);
+   RooChebychev bkg1("bkg1", "b1", xvar1, {a0});
+
+   RooRealVar c0("c0", "c0", -0.1, -1.0, 1.0);
+   RooChebychev bkg2("bkg2", "b2", xvar2, {c0});
+
+   RooRealVar s1("s1", "s1", 75.0, 0.0, 100000.0);
+   RooRealVar b1("b1", "b1", 25.0, 0.0, 100000.0);
+
+   // Extended models
+   RooAddPdf model1e("model1e", "", {gauss1, bkg1}, {s1, b1});
+
+   RooRealVar s2("s2", "s2", 50.0, 0.0, 100000.0);
+   RooRealVar b2("b2", "b2", 50.0, 0.0, 100000.0);
+
+   RooAddPdf model2e("model2e", "model2e", {gauss2, bkg2}, {s2, b2});
+
+   // Non-extended models
+   RooRealVar f1("f1", "f1", 0.75, 0.0, 1.0);
+   RooAddPdf model1n("model1n", "model1n", {gauss1, bkg1}, {f1});
+
+   RooRealVar f2("f2", "f2", 0.50, 0.0, 1.0);
+   RooAddPdf model2n("model2n", "model2n", RooArgList(gauss2, bkg2), RooArgList(f2));
+
+   // Case handling
+   enum class Case {
+      Gaussian,
+      NonExtended,
+      Extended
+   };
+
+   auto caseName = [](Case c) {
+      switch (c) {
+      case Case::Gaussian: return "Gaussian";
+      case Case::NonExtended: return "NonExtended";
+      case Case::Extended: return "Extended";
+      }
+      return "Unknown";
+   };
+
+   const std::vector<Case> cases = {Case::Gaussian, Case::NonExtended, Case::Extended};
+
+   // Helpers
+   auto integrateLastCurve = [](RooPlot *plot) {
+      const double xmin = plot->getPlotVar()->getMin();
+      const double xmax = plot->getPlotVar()->getMax();
+      return plot->getCurve()->average(xmin, xmax) * plot->getPlotVar()->numBins();
+   };
+
+   constexpr double tol = 0.01; // tolerate 1 % sampling error
+
+   // Test body
+   auto runCase = [&](Case c) {
+      SCOPED_TRACE(std::string("Case = ") + caseName(c));
+
+      RooAbsPdf *model1 = nullptr;
+      RooAbsPdf *model2 = nullptr;
+
+      switch (c) {
+      case Case::Gaussian:
+         model1 = &gauss1;
+         model2 = &gauss2;
+         break;
+      case Case::NonExtended:
+         model1 = &model1n;
+         model2 = &model2n;
+         break;
+      case Case::Extended:
+         model1 = &model1e;
+         model2 = &model2e;
+         break;
+      }
+
+      ASSERT_NE(model1, nullptr);
+      ASSERT_NE(model2, nullptr);
+
+      // Generate data
+      std::unique_ptr<RooDataSet> data1{model1->generate(xvar1, 10000)};
+      std::unique_ptr<RooDataSet> data2{model2->generate(xvar2, 1000)};
+
+      // Category
+      RooCategory sample("sample", "");
+      sample.defineType("Fit1");
+      sample.defineType("Fit2");
+
+      RooDataSet data("combinedData", "", {xvar1, xvar2}, Index(sample),
+                      Import({{"Fit1", data1.get()}, {"Fit2", data2.get()}}));
+
+      // Simultaneous PDF
+      RooSimultaneous sim_pdf("sim_pdf", "", sample);
+      sim_pdf.addPdf(*model1, "Fit1");
+      sim_pdf.addPdf(*model2, "Fit2");
+
+      std::unique_ptr<RooFitResult> result{sim_pdf.fitTo(data, Save(true), PrintLevel(-1))};
+
+      // Plot + checks
+      RooPlot *frame1 = xvar1.frame();
+      data.plotOn(frame1, Cut("sample==sample::Fit1"));
+      sim_pdf.plotOn(frame1, Slice(sample, "Fit1"), ProjWData(sample, data));
+
+      EXPECT_THAT(integrateLastCurve(frame1), RelativeNear(data1->sumEntries(), tol));
+
+      RooPlot *frame2 = xvar2.frame();
+      data.plotOn(frame2, Cut("sample==sample::Fit2"));
+      sim_pdf.plotOn(frame2, Slice(sample, "Fit2"), ProjWData(sample, data));
+
+      EXPECT_THAT(integrateLastCurve(frame2), RelativeNear(data2->sumEntries(), tol));
+   };
+
+   // Execute
+   for (Case c : cases) {
+      runCase(c);
+   }
+}
+
 /// JIRA ticket https://its.cern.ch/jira/browse/ROOT-7499
 /// Check that we can also generate Asimov datasets with non-integer weights
 /// via RooSimultaneous.
 TEST(RooSimultaneous, ExpectedDataWithNonIntegerWeights)
 {
+   RooHelpers::LocalChangeMsgLevel changeMsgLevel{RooFit::WARNING};
 
    RooWorkspace ws{"ws"};
    ws.factory("dummy_obs_a[0,1]");
diff --git a/tutorials/roofit/roofit/rf501_simultaneouspdf.C b/tutorials/roofit/roofit/rf501_simultaneouspdf.C
index a658fc66ea844..fd666c6b88039 100644
--- a/tutorials/roofit/roofit/rf501_simultaneouspdf.C
+++ b/tutorials/roofit/roofit/rf501_simultaneouspdf.C
@@ -102,29 +102,28 @@ void rf501_simultaneouspdf()
    RooPlot *frame1 = x.frame(Title("Physics sample"));
 
    // Plot all data tagged as physics sample
-   combData.plotOn(frame1, Cut("sample==sample::physics"));
+   std::unique_ptr<RooAbsData> slicedData1{combData.reduce(Cut("sample==sample::physics"))};
+   slicedData1->plotOn(frame1);
 
    // Plot "physics" slice of simultaneous pdf.
+   simPdf.getPdf("physics")->plotOn(frame1);
+   simPdf.getPdf("physics")->plotOn(frame1, Components("px"), LineStyle(kDashed));
+
+   // The same plot for the control sample slice. We do this with a different
+   // approach, using the data slice for the projection. This approach is more
+   // general, because you can plot sums of slices by using logical or in the
+   // Cut() command.
    // NBL You _must_ project the sample index category with data using ProjWData
    // as a RooSimultaneous makes no prediction on the shape in the index category
    // and can thus not be integrated.
    // In other words: Since the PDF doesn't know the number of events in the different
    // category states, it doesn't know how much of each component it has to project out.
    // This information is read from the data.
-   simPdf.plotOn(frame1, Slice(sample, "physics"), ProjWData(sample, combData));
-   simPdf.plotOn(frame1, Slice(sample, "physics"), Components("px"), ProjWData(sample, combData), LineStyle(kDashed));
-
-   // The same plot for the control sample slice. We do this with a different
-   // approach this time, for illustration purposes. Here, we are slicing the
-   // dataset and then use the data slice for the projection, because then the
-   // RooFit::Slice() becomes unnecessary. This approach is more general,
-   // because you can plot sums of slices by using logical or in the Cut()
-   // command.
    RooPlot *frame2 = x.frame(Bins(30), Title("Control sample"));
-   std::unique_ptr<RooAbsData> slicedData{combData.reduce(Cut("sample==sample::control"))};
-   slicedData->plotOn(frame2);
-   simPdf.plotOn(frame2, ProjWData(sample, *slicedData));
-   simPdf.plotOn(frame2, Components("px_ctl"), ProjWData(sample, *slicedData), LineStyle(kDashed));
+   std::unique_ptr<RooAbsData> slicedData2{combData.reduce(Cut("sample==sample::control"))};
+   slicedData2->plotOn(frame2);
+   simPdf.plotOn(frame2, ProjWData(sample, *slicedData2));
+   simPdf.plotOn(frame2, Components("px_ctl"), ProjWData(sample, *slicedData2), LineStyle(kDashed));
 
    // The same plot for all the phase space. Here, we can just use the original
    // combined dataset.
diff --git a/tutorials/roofit/roofit/rf501_simultaneouspdf.py b/tutorials/roofit/roofit/rf501_simultaneouspdf.py
index 4c645b91a922c..6f4f81da61807 100644
--- a/tutorials/roofit/roofit/rf501_simultaneouspdf.py
+++ b/tutorials/roofit/roofit/rf501_simultaneouspdf.py
@@ -13,7 +13,6 @@
 
 import ROOT
 
-
 # Create model for physics sample
 # -------------------------------------------------------------
 
@@ -96,28 +95,28 @@
 frame1 = x.frame(Title="Physics sample")
 
 # Plot all data tagged as physics sample
-combData.plotOn(frame1, Cut="sample==sample::physics")
+slicedData1 = combData.reduce(Cut="sample==sample::physics")
+slicedData1.plotOn(frame1)
 
 # Plot "physics" slice of simultaneous pdf.
-# NB: You *must* project the sample index category with data using ProjWData as
-# a RooSimultaneous makes no prediction on the shape in the index category and
-# can thus not be integrated. In other words: Since the PDF doesn't know the
-# number of events in the different category states, it doesn't know how much
-# of each component it has to project out. This info is read from the data.
-simPdf.plotOn(frame1, Slice=(sample, "physics"), ProjWData=(sample, combData))
-simPdf.plotOn(frame1, Slice=(sample, "physics"), Components="px", ProjWData=(sample, combData), LineStyle="--")
+simPdf.getPdf("physics").plotOn(frame1)
+simPdf.getPdf("physics").plotOn(frame1, Components="px", LineStyle="--")
 
 # The same plot for the control sample slice. We do this with a different
-# approach this time, for illustration purposes. Here, we are slicing the
-# dataset and then use the data slice for the projection, because then the
-# RooFit::Slice() becomes unnecessary. This approach is more general,
-# because you can plot sums of slices by using logical or in the Cut()
-# command.
+# approach, using the data slice for the projection. This approach is more
+# general, because you can plot sums of slices by using logical or in the
+# Cut() command.
+# NBL You _must_ project the sample index category with data using ProjWData
+# as a RooSimultaneous makes no prediction on the shape in the index category
+# and can thus not be integrated.
+# In other words: Since the PDF doesn't know the number of events in the different
+# category states, it doesn't know how much of each component it has to project out.
+# This information is read from the data.
 frame2 = x.frame(Title="Control sample")
-slicedData = combData.reduce(Cut="sample==sample::control")
-slicedData.plotOn(frame2)
-simPdf.plotOn(frame2, ProjWData=(sample, slicedData))
-simPdf.plotOn(frame2, Components="px_ctl", ProjWData=(sample, slicedData), LineStyle="--")
+slicedData2 = combData.reduce(Cut="sample==sample::control")
+slicedData2.plotOn(frame2)
+simPdf.plotOn(frame2, ProjWData=(sample, slicedData2))
+simPdf.plotOn(frame2, Components="px_ctl", ProjWData=(sample, slicedData2), LineStyle="--")
 
 # The same plot for all the phase space. Here, we can just use the original
 # combined dataset.