    /*****************************************************************************
     * Project: RooFit                                                           *
     * Package: RooFitCore                                                       *
     *    File: $Id: RooProdPdf.cc,v 1.60 2005/06/23 07:37:30 wverkerke Exp $
     * Authors:                                                                  *
     *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
     *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
     *                                                                           *
     * Copyright (c) 2000-2005, Regents of the University of California          *
     *                          and Stanford University. All rights reserved.    *
     *                                                                           *
     * Redistribution and use in source and binary forms,                        *
     * with or without modification, are permitted according to the terms        *
     * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
     *****************************************************************************/

    // -- CLASS DESCRIPTION [PDF] --
    // RooProdPdf is an efficient implementation of a product of PDFs of the form
    //
    //  PDF_1 * PDF_2 * ... * PDF_N
    //
    // PDFs may share observables. If that is the case any irreducable subset
    // of PDFS that share observables will be normalized with explicit numeric
    // integration as any built-in normalization will no longer be valid.
    //
    // Alternatively, products using conditional PDFs can be defined, e.g.
    //
    //    F(x|y) * G(y)
    //
    // meaning a pdf F(x) _given_ y and a PDF G(y). In this contruction F is only
    // normalized w.r.t x and G is normalized w.r.t y. The product in this construction
    // is properly normalized.
    //
    // If exactly one of the component PDFs supports extended likelihood fits, the
    // product will also be usable in extended mode, returning the number of expected
    // events from the extendable component PDF. The extendable component does not
    // have to appear in any specific place in the list.
    //

    #include "RooFit.h"

    #include "TIterator.h"
    #include "TIterator.h"
    #include "RooProdPdf.h"
    #include "RooRealProxy.h"
    #include "RooProdGenContext.h"
    #include "RooGenProdProj.h"
    #include "RooProduct.h"
    #include "RooNameReg.h"

    ClassImp(RooProdPdf)
    ;

     RooProdPdf::RooProdPdf(const char *name, const char *title, Double_t cutOff) :
      RooAbsPdf(name,title),
      _partListMgr(10),
      _partOwnedListMgr(10),
      _genCode(10),
      _cutOff(cutOff),
      _pdfList("_pdfList","List of PDFs",this),
      _extendedIndex(-1),
      _useDefaultGen(kFALSE)
    {
      // Dummy constructor
    }

     RooProdPdf::RooProdPdf(const char *name, const char *title,
    		       RooAbsPdf& pdf1, RooAbsPdf& pdf2, Double_t cutOff) :
      RooAbsPdf(name,title),
      _partListMgr(10),
      _partOwnedListMgr(10),
      _genCode(10),
      _cutOff(cutOff),
      _pdfList("_pdfList","List of PDFs",this),
      _pdfIter(_pdfList.createIterator()),
      _extendedIndex(-1),
      _useDefaultGen(kFALSE)
    {
      // Constructor with 2 PDFs (most frequent use case).
      //
      // The optional cutOff parameter can be used as a speed optimization if
      // one or more of the PDF have sizable regions with very small values,
      // which would pull the entire product of PDFs to zero in those regions.
      //
      // After each PDF multiplication, the running product is compared with
      // the cutOff parameter. If the running product is smaller than the
      // cutOff value, the product series is terminated and remaining PDFs
      // are not evaluated.
      //
      // There is no magic value of the cutOff, the user should experiment
      // to find the appropriate balance between speed and precision.
      // If a cutoff is specified, the PDFs most likely to be small should
      // be put first in the product. The default cutOff value is zero.

      _pdfList.add(pdf1) ;
      RooArgSet* nset1 = new RooArgSet("nset1") ;
      _pdfNSetList.Add(nset1) ;
      if (pdf1.canBeExtended()) {
        _extendedIndex = _pdfList.index(&pdf1) ;
      }

      _pdfList.add(pdf2) ;
      RooArgSet* nset2 = new RooArgSet("nset2") ;
      _pdfNSetList.Add(nset2) ;

      if (pdf2.canBeExtended()) {
        if (_extendedIndex>=0) {
          // Protect against multiple extended terms
          cout << "RooProdPdf::RooProdPdf(" << GetName()
    	   << ") multiple components with extended terms detected,"
    	   << " product will not be extendible." << endl ;
          _extendedIndex=-1 ;
        } else {
          _extendedIndex=_pdfList.index(&pdf2) ;
        }
      }
    }

     RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgList& pdfList, Double_t cutOff) :
      RooAbsPdf(name,title),
      _partListMgr(10),
      _partOwnedListMgr(10),
      _genCode(10),
      _cutOff(cutOff),
      _pdfList("_pdfList","List of PDFs",this),
      _pdfIter(_pdfList.createIterator()),
      _extendedIndex(-1),
      _useDefaultGen(kFALSE)
    {
      // Constructor from a list of PDFs
      //
      // The optional cutOff parameter can be used as a speed optimization if
      // one or more of the PDF have sizable regions with very small values,
      // which would pull the entire product of PDFs to zero in those regions.
      //
      // After each PDF multiplication, the running product is compared with
      // the cutOff parameter. If the running product is smaller than the
      // cutOff value, the product series is terminated and remaining PDFs
      // are not evaluated.
      //
      // There is no magic value of the cutOff, the user should experiment
      // to find the appropriate balance between speed and precision.
      // If a cutoff is specified, the PDFs most likely to be small should
      // be put first in the product. The default cutOff value is zero.

      TIterator* iter = pdfList.createIterator() ;
      RooAbsArg* arg ;
      Int_t numExtended(0) ;
      while((arg=(RooAbsArg*)iter->Next())) {
        RooAbsPdf* pdf = dynamic_cast<RooAbsPdf*>(arg) ;
        if (!pdf) {
          cout << "RooProdPdf::RooProdPdf(" << GetName() << ") list arg "
    	   << pdf->GetName() << " is not a PDF, ignored" << endl ;
          continue ;
        }
        _pdfList.add(*pdf) ;

        RooArgSet* nset = new RooArgSet("nset") ;
        _pdfNSetList.Add(nset) ;

        if (pdf->canBeExtended()) {
          _extendedIndex = _pdfList.index(pdf) ;
          numExtended++ ;
        }
      }

      // Protect against multiple extended terms
      if (numExtended>1) {
        cout << "RooProdPdf::RooProdPdf(" << GetName()
    	 << ") WARNING: multiple components with extended terms detected,"
    	 << " product will not be extendible." << endl ;
        _extendedIndex = -1 ;
      }

      delete iter ;
    }

     RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgSet& fullPdfSet,
    		       const RooCmdArg& arg1, const RooCmdArg& arg2,
    		       const RooCmdArg& arg3, const RooCmdArg& arg4,
    		       const RooCmdArg& arg5, const RooCmdArg& arg6,
    		       const RooCmdArg& arg7, const RooCmdArg& arg8) :
      RooAbsPdf(name,title),
      _partListMgr(10),
      _partOwnedListMgr(10),
      _genCode(10),
      _cutOff(0),
      _pdfList("_pdfList","List of PDFs",this),
      _pdfIter(_pdfList.createIterator()),
      _extendedIndex(-1),
      _useDefaultGen(kFALSE)
    {
      // Constructor from named argument list
      //
      // fullPdf -- Set of 'regular' PDFS that are normalized over all their observables
      // ConditionalPdf(pdfSet,depSet) -- Add PDF to product with condition that it
      //                                  only be normalized over specified observables
      //                                  any remaining observables will be conditional
      //                                  observables
      //
      //
      // For example, given a PDF F(x,y) and G(y)
      //
      // RooProdPdf("P","P",G,Partial(F,x)) will construct a 2-dimensional PDF as follows:
      //
      //   P(x,y) = G(y)/Int[y]G(y) * F(x,y)/Int[x]G(x,y)
      //
      // which is a well normalized and properly defined PDF, but different from the
      //
      //  P'(x,y) = F(x,y)*G(y) / Int[x,y] F(x,y)*G(y)
      //
      // In the former case the y distribution of P is identical to that of G, while
      // F only is used to determine the correlation between X and Y. In the latter
      // case the Y distribution is defined by the product of F and G.
      //
      // This P(x,y) construction is analoguous to generating events from F(x,y) with
      // a prototype dataset sampled from G(y)

      RooLinkedList l ;
      l.Add((TObject*)&arg1) ;  l.Add((TObject*)&arg2) ;
      l.Add((TObject*)&arg3) ;  l.Add((TObject*)&arg4) ;
      l.Add((TObject*)&arg5) ;  l.Add((TObject*)&arg6) ;
      l.Add((TObject*)&arg7) ;  l.Add((TObject*)&arg8) ;

      initializeFromCmdArgList(fullPdfSet,l) ;
    }

     RooProdPdf::RooProdPdf(const char* name, const char* title, const RooArgSet& fullPdfSet, const RooLinkedList& cmdArgList) :
      RooAbsPdf(name,title),
      _partListMgr(10),
      _partOwnedListMgr(10),
      _genCode(10),
      _cutOff(0),
      _pdfList("_pdfList","List of PDFs",this),
      _pdfIter(_pdfList.createIterator()),
      _extendedIndex(-1),
      _useDefaultGen(kFALSE)
    {
    //   cout << "RooProdPdf::ctor" << endl ;
    //   cmdArgList.Print("v") ;
      initializeFromCmdArgList(fullPdfSet, cmdArgList) ;
    }

     RooProdPdf::RooProdPdf(const RooProdPdf& other, const char* name) :
      RooAbsPdf(other,name),
      _partListMgr(other._partListMgr,kTRUE),
      _partOwnedListMgr(other._partOwnedListMgr,kTRUE),
      _genCode(other._genCode),
      _cutOff(other._cutOff),
      _pdfList("_pdfList",this,other._pdfList),
      _pdfIter(_pdfList.createIterator()),
      _extendedIndex(other._extendedIndex),
      _useDefaultGen(other._useDefaultGen)
    {
      // Copy constructor

      // Clone contents of normalizarion set list
      TIterator* iter = other._pdfNSetList.MakeIterator() ;
      RooArgSet* nset ;
      while((nset=(RooArgSet*)iter->Next())) {
        _pdfNSetList.Add(nset->snapshot()) ;
      }
      delete iter ;

    }

     void RooProdPdf::initializeFromCmdArgList(const RooArgSet& fullPdfSet, const RooLinkedList& l)
    {
      // Initialize RooProdPdf from a list of RooCmdArg configuration arguments

      // Process set of full PDFS
      TIterator* siter = fullPdfSet.createIterator() ;
      RooAbsPdf* pdf ;
      while((pdf=(RooAbsPdf*)siter->Next())) {
        _pdfList.add(*pdf) ;
        RooArgSet* nset1 = new RooArgSet("nset1") ;
        _pdfNSetList.Add(nset1) ;
      }
      delete siter ;

      // Process list of conditional PDFs
      TIterator* iter = l.MakeIterator() ;
      RooCmdArg* carg ;
      while((carg=(RooCmdArg*)iter->Next())) {
        if (!TString(carg->GetName()).CompareTo("Conditional")) {

          RooArgSet* pdfSet = (RooArgSet*) carg->getObject(0) ;
          RooArgSet* normSet = (RooArgSet*) carg->getObject(1) ;
          TIterator* siter = pdfSet->createIterator() ;
          RooAbsPdf* pdf ;
          while((pdf=(RooAbsPdf*)siter->Next())) {
    	_pdfList.add(*pdf) ;
    	_pdfNSetList.Add(normSet->snapshot()) ;
          }
          delete siter ;

        } else if (TString(carg->GetName()).CompareTo("")) {
          cout << "Unknown arg: " << carg->GetName() << endl ;
        }
      }
      delete iter ;
    }

     RooProdPdf::~RooProdPdf()
    {
      // Destructor
      Int_t i ;

      _pdfNSetList.Delete() ;
      for (i=0 ; i<10 ; i++) {
        _partNormListCache[i].Delete() ;
      }

      delete _pdfIter ;
    }

     Double_t RooProdPdf::getVal(const RooArgSet* set) const
    {
      _curNormSet = (RooArgSet*)set ;
      return RooAbsPdf::getVal(set) ;
    }

     Double_t RooProdPdf::evaluate() const
    {
      // Calculate current unnormalized value of object

      // WVE NSET may not exists, but getPartIntList may need it...

      Int_t code ;
      RooArgList *plist ;
      RooLinkedList *nlist ;
      getPartIntList(_curNormSet,0,plist,nlist,code) ;

      return calculate(plist,nlist) ;
    }

     Double_t RooProdPdf::calculate(const RooArgList* partIntList, const RooLinkedList* normSetList) const
    {
      // Calculate running product of pdfs, skipping factorized components
      RooAbsReal* partInt ;
      RooArgSet* normSet ;
      Double_t value(1.0) ;
      Int_t n = partIntList->getSize() ;

      Int_t i ;
      for (i=0 ; i<n ; i++) {
        partInt = ((RooAbsReal*)partIntList->at(i)) ;
        normSet = ((RooArgSet*)normSetList->At(i)) ;
        Double_t piVal = partInt->getVal(normSet->getSize()>0 ? normSet : 0) ;
        value *= piVal ;
    //     if (_verboseEval<0) {
    //       cout << "RPP:calc(" << GetName() << "): value *= " << piVal << " (" << partInt->GetName() << ") nset = " ; normSet->Print("1") ;
    //     }
        if (value<_cutOff) {
          //cout << "RooProdPdf::calculate(" << GetName() << ") calculation cut off after " << partInt->GetName() << endl ;
          break ;
        }
      }

      return value ;
    }

     void RooProdPdf::factorizeProduct(const RooArgSet& normSet, const RooArgSet& intSet,
    				  RooLinkedList& termList, RooLinkedList& normList,
    				  RooLinkedList& impDepList, RooLinkedList& crossDepList,
    				  RooLinkedList& intList) const
    {
      // Factorize product in irreducible terms for given choice of integration/normalization

    //   cout << "RooProdPdf::factorizeProduct(" << GetName() << ")" << endl ;
    //   cout << "   normSet  = " ; normSet.Print("1") ;
    //   cout << "    intSet  = " ; intSet.Print("1") ;

      _pdfIter->Reset() ;
      RooAbsPdf* pdf ;

      // List of all term dependents: normalization and imported
      RooLinkedList depAllList ;
      RooLinkedList depIntNoNormList ;

      // Setup lists for factorization terms and their dependents
      RooArgSet* term(0) ;
      RooArgSet* termNormDeps(0) ;
      RooArgSet* termAllDeps(0) ;
      RooArgSet* termIntDeps(0) ;
      RooArgSet* termIntNoNormDeps(0) ;
      TIterator* lIter = termList.MakeIterator() ;
      TIterator* ldIter = normList.MakeIterator() ;
      TIterator* laIter = depAllList.MakeIterator() ;
      TIterator* nIter = _pdfNSetList.MakeIterator() ;
      RooArgSet* pdfNSet ;

      // Loop over the PDFs
      while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
        pdfNSet = (RooArgSet*) nIter->Next() ;
        lIter->Reset() ;
        ldIter->Reset() ;
        laIter->Reset() ;

        // Reduce pdfNSet to actual dependents
        pdfNSet = pdf->getObservables(*pdfNSet) ;
    //     cout << "factorize: pdf = " << pdf->GetName() << " nset = " ; pdfNSet->Print("1") ;

        RooArgSet pdfNormDeps ; // Dependents to be normalized for the PDF
        RooArgSet pdfAllDeps ; // All dependents of this PDF ;

        // Make list of all dependents of this PDF
        RooArgSet* tmp = pdf->getObservables(normSet) ;
        pdfAllDeps.add(*tmp) ;
        delete tmp ;

        // Make list of normalization dependents for this PDF ;
        if (pdfNSet->getSize()>0) {
          RooArgSet* tmp = (RooArgSet*) pdfAllDeps.selectCommon(*pdfNSet) ;
          pdfNormDeps.add(*tmp) ;
          delete tmp ;
        } else {
          pdfNormDeps.add(pdfAllDeps) ;
        }

        RooArgSet* pdfIntSet = pdf->getObservables(intSet) ;

        RooArgSet pdfIntNoNormDeps(*pdfIntSet) ;
        pdfIntNoNormDeps.remove(pdfNormDeps,kTRUE,kTRUE) ;

    //     cout << "factorize: all pdf dependents = " ; pdfAllDeps.Print("1") ;
    //     cout << "factorize: normalized pdf dependents = " ; pdfNormDeps.Print("1") ;
    //     cout << "factorize: integrated pdf dependents = " ; pdfIntSet->Print("1") ;
    //     cout << "factorize: int/no-norm pdf dependents = " ; pdfIntNoNormDeps.Print("1") ;

        // Check if this PDF has dependents overlapping with one of the existing terms
        Bool_t done(kFALSE) ;
        while((term=(RooArgSet*)lIter->Next())) {
          termNormDeps=(RooArgSet*)ldIter->Next() ;
          termAllDeps=(RooArgSet*)laIter->Next() ;

          // PDF should be added to existing term if
          // 1) It has overlapping normalization dependents with any other PDF in existing term
          // 2) It has overlapping dependents of any class for which integration is requested

          Bool_t normOverlap = pdfNormDeps.overlaps(*termNormDeps)  ;
          // Bool_t intOverlap =  pdfIntSet->overlaps(*termAllDeps) ;

          //if (normOverlap || intOverlap) {
          if (normOverlap) {
    // 	cout << "adding pdf " << pdf->GetName() << " to term " ; term->Print("1") ;
    // 	cout << "   with termNormDeps " ; termNormDeps->Print("1") ;
    	term->add(*pdf) ;
    	termNormDeps->add(pdfNormDeps,kFALSE) ;
    	termAllDeps->add(pdfAllDeps,kFALSE) ;
    	termIntDeps->add(*pdfIntSet,kFALSE) ;
    	termIntNoNormDeps->add(pdfIntNoNormDeps,kFALSE) ;
    	done = kTRUE ;
          }
        }

        // If not, create a new term
        if (!done) {
    //       cout << "creating new term" << endl ;
          term = new RooArgSet("term") ;
          termNormDeps = new RooArgSet("termNormDeps") ;
          termAllDeps = new RooArgSet("termAllDeps") ;
          termIntDeps = new RooArgSet("termIntDeps") ;
          termIntNoNormDeps = new RooArgSet("termIntNoNormDeps") ;

          term->add(*pdf) ;
          termNormDeps->add(pdfNormDeps,kFALSE) ;
          termAllDeps->add(pdfAllDeps,kFALSE) ;
          termIntDeps->add(*pdfIntSet,kFALSE) ;
          termIntNoNormDeps->add(pdfIntNoNormDeps,kFALSE) ;

          termList.Add(term) ;
          normList.Add(termNormDeps) ;
          depAllList.Add(termAllDeps) ;
          intList.Add(termIntDeps) ;
          depIntNoNormList.Add(termIntNoNormDeps) ;
        }

        // We own the reduced version of pdfNSet
        delete pdfNSet ;
        delete pdfIntSet ;
      }

      // Loop over list of terms again to determine 'imported' observables
      lIter->Reset() ;
      ldIter->Reset() ;
      laIter->Reset() ;
      TIterator* innIter = depIntNoNormList.MakeIterator() ;

      while((term=(RooArgSet*)lIter->Next())) {
        RooArgSet* normDeps = (RooArgSet*) ldIter->Next() ;
        RooArgSet* allDeps = (RooArgSet*) laIter->Next() ;
        RooArgSet* intNoNormDeps = (RooArgSet*) innIter->Next() ;

        // Make list of wholly imported dependents
        RooArgSet impDeps(*allDeps) ;
        impDeps.remove(*normDeps,kTRUE,kTRUE) ;
        impDepList.Add(impDeps.snapshot()) ;

        // Make list of cross dependents (term is self contained for these dependents,
        // but components import dependents from other components)
        RooArgSet* crossDeps = (RooArgSet*) intNoNormDeps->selectCommon(*normDeps) ;
        crossDepList.Add(crossDeps->snapshot()) ;
        delete crossDeps ;
      }

    //    lIter->Reset() ;
    //    cout << "list of terms:" << endl ;
    //    while((term=(RooArgSet*)lIter->Next())) {
    //      term->Print("1") ;
    //    }

      depAllList.Delete() ;
      depIntNoNormList.Delete() ;

      delete nIter ;
      delete lIter ;
      delete ldIter ;
      delete laIter ;
      delete innIter ;

      return ;
    }

     void RooProdPdf::getPartIntList(const RooArgSet* nset, const RooArgSet* iset,
    				pRooArgList& partList, pRooLinkedList& nsetList,
    				Int_t& code, const char* isetRangeName) const
    {
      // Check if this configuration was created before
      Int_t sterileIdx(-1) ;
      RooArgList* partIntList = _partListMgr.getNormList(this,nset,iset,&sterileIdx,RooNameReg::ptr(isetRangeName)) ;
      if (partIntList) {
        code = _partListMgr.lastIndex() ;

        partList = partIntList ;
        nsetList = _partNormListCache+code ;

        return ;
      }

      // Create containers for partial integral components to be generated
      partIntList = new RooArgList("partIntList") ;
      RooArgList* partIntOwnedList = new RooArgList("partIntList") ;
      RooLinkedList* partIntNormList = new RooLinkedList ;

      // Factorize the product in irreducible terms for this nset
      RooLinkedList terms, norms, imp, ints, cross ;
      factorizeProduct(nset?(*nset):RooArgSet(),iset?(*iset):RooArgSet(),terms,norms,imp,cross,ints) ;

      RooArgSet *norm, *integ, *xdeps ;

      // Group irriducible terms that need to be (partially) integrated together
      RooLinkedList groupedList ;
      RooArgSet outerIntDeps ;
      groupProductTerms(groupedList,outerIntDeps,terms,norms,imp,ints,cross) ;
      TIterator* gIter = groupedList.MakeIterator() ;
      RooLinkedList* group ;

      while((group=(RooLinkedList*)gIter->Next())) {

    //     group->Print("1") ;

        if (group->GetSize()==1) {
          RooArgSet* term = (RooArgSet*) group->At(0) ;

          Int_t termIdx = terms.IndexOf(term) ;
          norm=(RooArgSet*) norms.At(termIdx) ;
          integ=(RooArgSet*) ints.At(termIdx) ;
          xdeps=(RooArgSet*)cross.At(termIdx) ;

          RooArgSet termNSet, termISet, termXSet ;

          // Take list of normalization, integrated dependents from factorization algorithm
          termISet.add(*integ) ;
          termNSet.add(*norm) ;

          // Cross-imported integrated dependents
          termXSet.add(*xdeps) ;

          // Add prefab term to partIntList.
          Bool_t isOwned ;
          RooAbsReal* func = processProductTerm(nset,iset,isetRangeName,term,termNSet,termISet,isOwned) ;
          if (func) {
    	partIntList->add(*func) ;
    	if (isOwned) partIntOwnedList->addOwned(*func) ;
    	partIntNormList->Add(norm->snapshot()) ;
          }
        } else {

    //       cout << "prod: composite group encountered!" << endl ;

          RooArgSet compTermSet, compTermNorm ;
          TIterator* tIter = group->MakeIterator() ;
          RooArgSet* term ;
          while((term=(RooArgSet*)tIter->Next())) {

    	Int_t termIdx = terms.IndexOf(term) ;
    	norm=(RooArgSet*) norms.At(termIdx) ;
    	integ=(RooArgSet*) ints.At(termIdx) ;
    	xdeps=(RooArgSet*)cross.At(termIdx) ;

    	RooArgSet termNSet, termISet, termXSet ;
    	termISet.add(*integ) ;
    	termNSet.add(*norm) ;
    	termXSet.add(*xdeps) ;

    	// Remove outer integration dependents from termISet
    	termISet.remove(outerIntDeps,kTRUE,kTRUE) ;
    // 	cout << "termISet = " ; termISet.Print("1") ;

    	Bool_t isOwned ;
    	RooAbsReal* func = processProductTerm(nset,iset,isetRangeName,term,termNSet,termISet,isOwned,kTRUE) ;
    // 	cout << "created composite term component " << func->GetName() << endl ;
    	if (func) {
    	  compTermSet.add(*func) ;
    	  if (isOwned) partIntOwnedList->addOwned(*func) ;
    	  compTermNorm.add(*norm,kFALSE) ;
    	}
          }

    //       cout << "constructing special composite product" << endl ;
    //       cout << "compTermSet = " ; compTermSet.Print("1") ;

          // WVE this doesn't work (yet) for some reason
          //const char* name = makeRGPPName("SPECPROJ_",compTermSet,outerIntDeps,RooArgSet()) ;
          //RooAbsReal* func = new RooGenProdProj(name,name,compTermSet,outerIntDeps,RooArgSet()) ;
          //partIntList->add(*func) ;
          //partIntOwnedList->addOwned(*func) ;

          // WVE but this does, so we'll keep it for the moment
          const char* prodname = makeRGPPName("SPECPROD_",compTermSet,outerIntDeps,RooArgSet(),isetRangeName) ;
          RooProduct* prodtmp = new RooProduct(prodname,prodname,compTermSet) ;
          partIntOwnedList->addOwned(*prodtmp) ;

          const char* intname = makeRGPPName("SPECINT_",compTermSet,outerIntDeps,RooArgSet(),isetRangeName) ;
          RooRealIntegral* inttmp = new RooRealIntegral(intname,intname,*prodtmp,outerIntDeps) ;
          partIntOwnedList->addOwned(*inttmp) ;
          partIntList->add(*inttmp);

          partIntNormList->Add(compTermNorm.snapshot()) ;

          delete tIter ;
        }
      }
      delete gIter ;

      // Store the partial integral list and return the assigned code ;
      code = _partListMgr.setNormList(this,nset,iset,partIntList,RooNameReg::ptr(isetRangeName)) ;
      _partOwnedListMgr.setNormList(this,nset,iset,partIntOwnedList,RooNameReg::ptr(isetRangeName)) ;

      // Store the normalization set list in the cache using the index code of _partListMgr
      _partNormListCache[code].Delete() ;
      // We own content of 'norms' and transfer ownership to the cache
      _partNormListCache[code] = *partIntNormList ;

      // WVE PRINT SECTION
    //       cout << "RooProdPdf::getPIL(" << GetName() << "," << this << ") creating new configuration with code " << code << endl
    //            << "    nset = " ; if (nset) nset->Print("1") ; else cout << "<none>" << endl ;
    //       cout << "    iset" << (isetRangeName?(Form("[%s]",isetRangeName)):"") << " = " ; if (iset) iset->Print("1") ; else cout << "<none>" << endl ;
    //       cout << "    Partial Integral List:" ;
    //       partIntList->Print("1") ;
    //       cout << "    Partial Owned Integral List:" ;
    //       partIntOwnedList->Print("1") ;
    //       cout << "    Normalization list: " << endl ; partIntNormList->Print("1") ;
    //       cout << endl  ;

      partList =  partIntList ;
      nsetList = _partNormListCache+code ;

      // We own contents of all lists filled by factorizeProduct()
      terms.Delete() ;
      ints.Delete() ;
      imp.Delete() ;
      norms.Delete() ;
      cross.Delete() ;
    }

     void RooProdPdf::groupProductTerms(RooLinkedList& groupedTerms, RooArgSet& outerIntDeps,
    				   const RooLinkedList& terms, const RooLinkedList& norms,
    				   const RooLinkedList& imps, const RooLinkedList& ints, const RooLinkedList& /*cross*/) const
    {
      // Start out with each term in its own group
      TIterator* tIter = terms.MakeIterator() ;
      RooArgSet* term ;
      while((term=(RooArgSet*)tIter->Next())) {
        RooLinkedList* group = new RooLinkedList ;
        group->Add(term) ;
        groupedTerms.Add(group) ;
      }
      delete tIter ;

      // Make list of imported dependents that occur in any term
      RooArgSet allImpDeps ;
      TIterator* iIter = imps.MakeIterator() ;
      RooArgSet *impDeps ;
      while((impDeps=(RooArgSet*)iIter->Next())) {
        allImpDeps.add(*impDeps,kFALSE) ;
      }
      delete iIter ;

      // Make list of integrated dependents that occur in any term
      RooArgSet allIntDeps ;
      iIter = ints.MakeIterator() ;
      RooArgSet *intDeps ;
      while((intDeps=(RooArgSet*)iIter->Next())) {
        allIntDeps.add(*intDeps,kFALSE) ;
      }
      delete iIter ;

    //   cout << "Complete lists of imported dependens" ; allImpDeps.Print("1") ;
    //   cout << "Complete lists of integrated dependens" ; allIntDeps.Print("1") ;

      RooArgSet* tmp = (RooArgSet*) allIntDeps.selectCommon(allImpDeps) ;
      outerIntDeps.removeAll() ;
      outerIntDeps.add(*tmp) ;
      delete tmp ;
    //   cout << "Integrated Dependents that need 'outer' integration: " ; outerIntDeps.Print("1") ;

      // Now iteratively merge groups that should be (partially) integrated together
      TIterator* oidIter = outerIntDeps.createIterator() ;
      TIterator* glIter = groupedTerms.MakeIterator() ;
      RooAbsArg* outerIntDep ;
      while ((outerIntDep =(RooAbsArg*)oidIter->Next())) {

    //     cout << "merge for outerIntDep " << outerIntDep->GetName() << endl ;

        // Collect groups that feature this dependent
        RooLinkedList* newGroup = 0 ;

        // Loop over groups
        RooLinkedList* group ;
        glIter->Reset() ;
        Bool_t needMerge = kFALSE ;
        while((group=(RooLinkedList*)glIter->Next())) {

    //       cout << "considering the following group:" << endl ;
    //       group->Print("1") ;
    //       cout << "----" << endl ;

          // See if any term in this group depends in any ay on outerDepInt
          RooArgSet* term ;
          TIterator* tIter = group->MakeIterator() ;
          while((term=(RooArgSet*)tIter->Next())) {

    	Int_t termIdx = terms.IndexOf(term) ;
    	RooArgSet* termNormDeps = (RooArgSet*) norms.At(termIdx) ;
    	RooArgSet* termIntDeps = (RooArgSet*) ints.At(termIdx) ;
    	RooArgSet* termImpDeps = (RooArgSet*) imps.At(termIdx) ;

    	if (termNormDeps->contains(*outerIntDep) ||
    	    termIntDeps->contains(*outerIntDep) ||
    	    termImpDeps->contains(*outerIntDep)) {
    	  needMerge = kTRUE ;
    	}

          }

          if (needMerge) {
    // 	cout << "group needs to be merged" << endl ;
    	// Create composite group if not yet existing
    	if (newGroup==0) {
    	  newGroup = new RooLinkedList ;
    	}

    	// Add terms of this group to new term
    	tIter->Reset() ;
    	while((term=(RooArgSet*)tIter->Next())) {
    // 	  cout << "transferring group term to new merged group" << endl ;
    	  newGroup->Add(term) ;
    	}

    	// Remove this group from list and delete it (but not its contents)
    // 	cout << "removing group from groupList" << endl ;
    	groupedTerms.Remove(group) ;
    	delete group ;
          }
          delete tIter ;
        }
        // If a new group has been created to merge terms dependent on current outerIntDep, add it to group list
        if (newGroup) {
    //       cout << "adding merged group to group list" << endl ;
          groupedTerms.Add(newGroup) ;
    //       newGroup->Print("1") ;
    //       cout << "..." << endl ;
        }

      }

      delete glIter ;
      delete oidIter ;
    }

     RooAbsReal* RooProdPdf::processProductTerm(const RooArgSet* nset, const RooArgSet* iset, const char* isetRangeName,
    				           const RooArgSet* term,const RooArgSet& termNSet, const RooArgSet& termISet,
    				           Bool_t& isOwned, Bool_t forceWrap) const
    {
      // CASE I: factorizing term: term is integrated over all normalizing observables
      // -----------------------------------------------------------------------------
      // Check if all observbales of this term are integrated. If so the term cancels
      if (termNSet.getSize()>0 && termNSet.getSize()==termISet.getSize() && isetRangeName==0) {

        // Term factorizes
        return 0 ;
      }

      // CASE II: Dropped terms: if term is entirely unnormalized, it should be dropped
      // ------------------------------------------------------------------------------
      if (nset && termNSet.getSize()==0) {

        // Drop terms that are not asked to be normalized
        return 0 ;
      }

      if (iset && termISet.getSize()>0) {
        if (term->getSize()==1) {

          // CASE IIIa: Normalized and partially integrated single PDF term
          //---------------------------------------------------------------

          TIterator* pIter = term->createIterator() ;
          RooAbsPdf* pdf = (RooAbsPdf*) pIter->Next() ;
          delete pIter ;

          RooAbsReal* partInt = pdf->createIntegral(termISet,termNSet,isetRangeName) ;
          partInt->setOperMode(operMode()) ;

          isOwned=kTRUE ;
          return partInt ;

        } else {

          // CASE IIIb: Normalized and partially integrated composite PDF term
          //---------------------------------------------------------------

          // Use auxiliary class RooGenProdProj to calculate this term
          const char* name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
          RooAbsReal* partInt = new RooGenProdProj(name,name,*term,termISet,termNSet,isetRangeName) ;
          partInt->setOperMode(operMode()) ;

          isOwned=kTRUE ;
          return partInt ;
        }
      }

      // CASE IVa: Normalized non-integrated composite PDF term
      // -------------------------------------------------------
      if (nset && nset->getSize()>0 && term->getSize()>1) {
        // Composite term needs normalized integration
        const char* name = makeRGPPName("GENPROJ_",*term,termISet,termNSet,isetRangeName) ;
        RooAbsReal* partInt = new RooGenProdProj(name,name,*term,termISet,termNSet,isetRangeName) ;
        partInt->setOperMode(operMode()) ;

        isOwned=kTRUE ;
        return partInt ;
      }

      // CASE IVb: Normalized, non-integrated single PDF term
      // -----------------------------------------------------
      TIterator* pIter = term->createIterator() ;
      RooAbsPdf* pdf ;
      while((pdf=(RooAbsPdf*)pIter->Next())) {

        if (forceWrap) {

          // Construct representative name of normalization wrapper
          TString name(pdf->GetName()) ;
          name.Append("_NORM[") ;
          TIterator* nIter = termNSet.createIterator() ;
          RooAbsArg* arg ;
          Bool_t first(kTRUE) ;
          while((arg=(RooAbsArg*)nIter->Next())) {
    	if (!first) {
    	  name.Append(",") ;
    	} else {
    	  first=kFALSE ;
    	}
    	name.Append(arg->GetName()) ;
          }
          name.Append("]") ;
          delete nIter ;

          RooAbsReal* partInt = new RooRealIntegral(name.Data(),name.Data(),*pdf,RooArgSet(),&termNSet) ;
          isOwned=kTRUE ;
          return partInt ;

        } else {
          isOwned=kFALSE ;
          return pdf  ;
        }
      }
      delete pIter ;

      cout << "RooProdPdf::processProductTerm(" << GetName() << ") unidentified term!!!" << endl ;
      return 0;
    }

     const char* RooProdPdf::makeRGPPName(const char* pfx, const RooArgSet& term, const RooArgSet& iset, const RooArgSet& nset, const char* isetRangeName) const
    {
      // Make an appropriate name for a RooGenProdProj object in getPartIntList()

      static TString pname ;
      pname = pfx ;

      TIterator* pIter = term.createIterator() ;

      // Encode component names
      Bool_t first(kTRUE) ;
      RooAbsPdf* pdf ;
      while((pdf=(RooAbsPdf*)pIter->Next())) {
        if (first) {
          first = kFALSE ;
        } else {
          pname.Append("_X_") ;
        }
        pname.Append(pdf->GetName()) ;
      }
      delete pIter ;

      pname.Append(integralNameSuffix(iset,&nset,isetRangeName)) ;

      return pname.Data() ;
    }

     Bool_t RooProdPdf::forceAnalyticalInt(const RooAbsArg& /*dep*/) const
    {
      return kTRUE ;
    }

     Int_t RooProdPdf::getAnalyticalIntegralWN(RooArgSet& allVars, RooArgSet& analVars,
    					  const RooArgSet* normSet, const char* rangeName) const
    {
      // Determine which part (if any) of given integral can be performed analytically.
      // If any analytical integration is possible, return integration scenario code.
      //
      // RooProdPdf implements two strategies in implementing analytical integrals
      //
      // First, PDF components whose entire set of dependents are requested to be integrated
      // can be dropped from the product, as they will integrate out to 1 by construction
      //:
      // Second, RooProdPdf queries each remaining component PDF for its analytical integration
      // capability of the requested set ('allVars'). It finds the largest common set of variables
      // that can be integrated by all remaining components. If such a set exists, it reconfirms that
      // each component is capable of analytically integrating the common set, and combines the components
      // individual integration codes into a single integration code valid for RooProdPdf.

      if (_forceNumInt) return 0 ;

      // Declare that we can analytically integrate all requested observables
      analVars.add(allVars) ;

      // Retrieve (or create) the required partial integral list
      Int_t code ;
      RooArgList *plist ;
      RooLinkedList *nlist ;
      getPartIntList(normSet,&allVars,plist,nlist,code,rangeName) ;

      return code+1 ;
    }

     Double_t RooProdPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const
    {
      // Return analytical integral defined by given scenario code

      // No integration scenario
      if (code==0) {
        return getVal(normSet) ;
      }

      // WVE needs adaptation for rangename feature

      // Partial integration scenarios
      RooArgList* partIntList = _partListMgr.getNormListByIndex(code-1) ;

      // If cache has been sterilized, revive this slot
      if (partIntList==0) {
        RooArgSet* vars = getParameters(RooArgSet()) ;
        RooArgSet* nset = _partListMgr.nameSet1ByIndex(code-1)->select(*vars) ;
        RooArgSet* iset = _partListMgr.nameSet2ByIndex(code-1)->select(*vars) ;

        Int_t code2(-1) ;
        RooLinkedList *nlist ;
        getPartIntList(nset,iset,partIntList,nlist,code2,rangeName) ;

        delete vars ;
        delete nset ;
        delete iset ;
      }
      Double_t val = calculate(partIntList,_partNormListCache+code-1) ;

      //cout << "RPP::aIWN(" << GetName() << ") value = " << val << endl ;
      return val ;
    }

     Bool_t RooProdPdf::checkObservables(const RooArgSet* /*nset*/) const
    {
      // Check that none of the PDFs have overlapping dependents
      return kFALSE ;

    //   Bool_t ret(kFALSE) ;

    //   _pdfIter->Reset() ;
    //   RooAbsPdf* pdf, *pdf2 ;
    //   TIterator* iter2 = _pdfList.createIterator() ;
    //   while((pdf = (RooAbsPdf*)_pdfIter->Next())) {
    //     *iter2 = *_pdfIter ;
    //     while((pdf2 = (RooAbsPdf*)iter2->Next())) {
    //       if (pdf->observableOverlaps(nset,*pdf2)) {
    // 	cout << "RooProdPdf::checkObservables(" << GetName() << "): ERROR: PDFs " << pdf->GetName()
    // 	     << " and " << pdf2->GetName() << " have one or more dependents in common" << endl ;
    // 	ret = kTRUE ;
    //       }
    //     }
    //   }
    //   delete iter2 ;
    //   return ret ;
    }

     RooAbsPdf::ExtendMode RooProdPdf::extendMode() const
    {
      return (_extendedIndex>=0) ? ((RooAbsPdf*)_pdfList.at(_extendedIndex))->extendMode() : CanNotBeExtended ;
    }

     Double_t RooProdPdf::expectedEvents(const RooArgSet* nset) const
    {
      assert(_extendedIndex>=0) ;
      return ((RooAbsPdf*)_pdfList.at(_extendedIndex))->expectedEvents(nset) ;
    }

     RooAbsGenContext* RooProdPdf::genContext(const RooArgSet &vars, const RooDataSet *prototype,
    					 const RooArgSet* auxProto, Bool_t verbose) const
    {
      if (_useDefaultGen) return RooAbsPdf::genContext(vars,prototype,auxProto,verbose) ;
      return new RooProdGenContext(*this,vars,prototype,auxProto,verbose) ;
    }

     Int_t RooProdPdf::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t staticInitOK) const
    {
      if (!_useDefaultGen) return 0 ;

      // Find the subset directVars that only depend on a single PDF in the product
      RooArgSet directSafe ;
      TIterator* dIter = directVars.createIterator() ;
      RooAbsArg* arg ;
      while((arg=(RooAbsArg*)dIter->Next())) {
        if (isDirectGenSafe(*arg)) directSafe.add(*arg) ;
      }
      delete dIter ;

      // Now find direct integrator for relevant components ;
      _pdfIter->Reset() ;
      RooAbsPdf* pdf ;
      Int_t code[64], n(0) ;
      while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
        RooArgSet pdfDirect ;
        code[n] = pdf->getGenerator(directSafe,pdfDirect,staticInitOK) ;
        if (code[n]!=0) {
          generateVars.add(pdfDirect) ;
        }
        n++ ;
      }

      if (generateVars.getSize()>0) {
        Int_t masterCode = _genCode.store(code,n) ;
        return masterCode+1 ;
      } else {
        return 0 ;
      }
    }

     void RooProdPdf::initGenerator(Int_t code)
    {
      if (!_useDefaultGen) return ;

      const Int_t* codeList = _genCode.retrieve(code-1) ;
      _pdfIter->Reset() ;
      RooAbsPdf* pdf ;
      Int_t i(0) ;
      while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
        if (codeList[i]!=0) {
          pdf->initGenerator(codeList[i]) ;
        }
        i++ ;
      }
    }

     void RooProdPdf::generateEvent(Int_t code)
    {
      if (!_useDefaultGen) return ;

      const Int_t* codeList = _genCode.retrieve(code-1) ;
      _pdfIter->Reset() ;
      RooAbsPdf* pdf ;
      Int_t i(0) ;
      while((pdf=(RooAbsPdf*)_pdfIter->Next())) {
        if (codeList[i]!=0) {
          pdf->generateEvent(codeList[i]) ;
        }
        i++ ;
      }

    }

     void RooProdPdf::operModeHook()
    {
      Int_t i ;
      for (i=0 ; i<_partListMgr.cacheSize() ; i++) {
        RooArgList* plist = _partOwnedListMgr.getNormListByIndex(i) ;
        if (plist) {
          TIterator* iter = plist->createIterator() ;
          RooAbsArg* arg ;
          while((arg=(RooAbsArg*)iter->Next())) {
    	arg->setOperMode(_operMode) ;
          }
          delete iter ;
        }
      }
      return ;
    }

     void RooProdPdf::clearCache()
    {
    //   cout << "RooProdPdf::clearCache(" << GetName() << ") sterilizing cache content" << endl ;
      Int_t i ;
      for (i=0 ; i<10 ; i++) {
        _partNormListCache[i].Delete() ;
      }
      _partListMgr.sterilize()  ;
      _partOwnedListMgr.sterilize()  ;
    }

     Bool_t RooProdPdf::redirectServersHook(const RooAbsCollection& /*newServerList*/, Bool_t /*mustReplaceAll*/,
    				       Bool_t /*nameChange*/, Bool_t /*isRecursive*/)
    {
      // Throw away cache, as figuring out redirections on the cache is an unsolvable problem.
      clearCache() ;
      return kFALSE ;
    }

     void RooProdPdf::printCompactTreeHook(ostream& os, const char* indent)
    {
      Int_t i ;
      os << indent << "RooProdPdf begin partial integral cache" << endl ;

      for (i=0 ; i<_partListMgr.cacheSize() ; i++) {
        RooArgList* plist = _partListMgr.getNormListByIndex(i) ;
        if (plist) {
          TIterator* iter = plist->createIterator() ;
          RooAbsArg* arg ;
          TString indent2(indent) ;
          indent2 += Form("[%d] ",i) ;
          while((arg=(RooAbsArg*)iter->Next())) {
    	arg->printCompactTree(os,indent2) ;
          }
          delete iter ;
        }
      }

      os << indent << "RooProdPdf end partial integral cache" << endl ;
    }

     Bool_t RooProdPdf::isDirectGenSafe(const RooAbsArg& arg) const
    {
      // Only override base class behaviour if default generator method is enabled
      if (!_useDefaultGen) return RooAbsPdf::isDirectGenSafe(arg) ;

      // Argument may appear in only one PDF component
      _pdfIter->Reset() ;
      RooAbsPdf* pdf, *thePdf(0) ;
      while((pdf=(RooAbsPdf*)_pdfIter->Next())) {

        if (pdf->dependsOn(arg)) {
          // Found PDF depending on arg

          // If multiple PDFs depend on arg directGen is not safe
          if (thePdf) return kFALSE ;

          thePdf = pdf ;
        }
      }
      // Forward call to relevant component PDF
      return thePdf?(thePdf->isDirectGenSafe(arg)):kFALSE ;
    }

     RooArgSet* RooProdPdf::findPdfNSet(RooAbsPdf& pdf) const
    {
      // Look up user specified normalization set for given input PDF component

      Int_t idx = _pdfList.index(&pdf) ;
      if (idx<0) return 0 ;
      return (RooArgSet*) _pdfNSetList.At(idx) ;
    }