ROOT logo
// @(#)root/roostats:$Id: SPlot.cxx 29107 2009-06-19 14:26:42Z moneta $
// Author: Kyle Cranmer   28/07/2008

/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

/*****************************************************************************
 * Project: RooStats
 * Package: RooFit/RooStats  
 *
 * Authors:                     
 *   Original code from M. Pivk as part of MLFit package from BaBar.
 *   Modifications:
 *     Giacinto Piacquadio, Maurizio Pierini: modifications for new RooFit version
 *     George H. Lewis, Kyle Cranmer: generalized for weighted events
 * 
 * Porting to RooStats (with permission) by Kyle Cranmer, July 2008
 *
 *****************************************************************************/



//_________________________________________________
//BEGIN_HTML
// This class calculates sWeights used to create an sPlot.  
// The code is based on 
// ``SPlot: A statistical tool to unfold data distributions,'' 
//  Nucl. Instrum. Meth. A 555, 356 (2005) 
//  [arXiv:physics/0402083].
//
// An SPlot gives us  the distribution of some variable, x in our 
// data sample for a given species (eg. signal or background).  
// The result is similar to a likelihood projection plot, but no cuts are made, 
//  so every event contributes to the distribution.
//
// [Usage]
// To use this class, you first must have a pdf that includes
// yields for (possibly several) different species.
// Create an instance of the class by supplying a data set,
// the pdf, and a list of the yield variables.  The SPlot Class
// will calculate SWeights and include these as columns in the RooDataSet.
//END_HTML
//

#include <vector>
#include <map>

#include "RooStats/SPlot.h"
#include "RooAbsPdf.h"
#include "RooDataSet.h"
#include "RooRealVar.h"
#include "RooGlobalFunc.h"
#include "TTree.h"
#include "RooStats/RooStatsUtils.h" 


#include "TMatrixD.h"


ClassImp(RooStats::SPlot) ;

using namespace RooStats;



//__________________________________________________________________
SPlot::~SPlot()
{
  if(fSData)
    delete fSData;

}

//____________________________________________________________________
SPlot::SPlot():
  TNamed()
{
  // Default constructor

  RooArgList Args;

  fSWeightVars = Args;

  fSData = NULL;

}

//_________________________________________________________________
SPlot::SPlot(const char* name, const char* title):
  TNamed(name, title)
{

  RooArgList Args;

  fSWeightVars = Args;

  fSData = NULL;

}

//___________________________________________________________________
SPlot::SPlot(const char* name, const char* title, const RooDataSet &data):
  TNamed(name, title)
{
  //Constructor from a RooDataSet
  //No sWeighted variables are present
  
  RooArgList Args;
  
  fSWeightVars = Args;

  fSData = (RooDataSet*) &data;
}


//____________________________________________________________________
SPlot::SPlot(const SPlot &other):
  TNamed(other)
{
  // Copy Constructor from another SPlot
  
  RooArgList Args = (RooArgList) other.GetSWeightVars();
  
  fSWeightVars.addClone(Args);

  fSData = (RooDataSet*) other.GetSDataSet();
  
}


//______________________________________________________________________
SPlot::SPlot(const char* name, const char* title, RooDataSet& data, RooAbsPdf* pdf, 
	     const RooArgList &yieldsList, const RooArgSet &projDeps, 
	     bool includeWeights, bool cloneData, const char* newName):
  TNamed(name, title)
{
  if(cloneData == 1)
    fSData = (RooDataSet*) data.Clone(newName);
  else
    fSData = (RooDataSet*) &data;

  //Construct a new SPlot class,
  //calculate sWeights, and include them
  //in the RooDataSet of this class.

  this->AddSWeight(pdf, yieldsList, projDeps, includeWeights);
}


//___________________________________________________________________________
RooDataSet* SPlot::SetSData(RooDataSet* data)
{
  if(data)    {
    fSData = (RooDataSet*) data;
    return fSData;
  }  else
    return NULL;
}

//__________________________________________________________________________
RooDataSet* SPlot::GetSDataSet() const
{
  return fSData;
}  

//____________________________________________________________________________
Double_t SPlot::GetSWeight(Int_t numEvent, const char* sVariable) const
{

  if(numEvent > fSData->numEntries() )
    {
      coutE(InputArguments)  << "Invalid Entry Number" << endl;
      return -1;
    }

  if(numEvent < 0)
    {
      coutE(InputArguments)  << "Invalid Entry Number" << endl;
      return -1;    
    }

  Double_t totalYield = 0;

  std::string varname(sVariable);
  varname += "_sw";


  if(fSWeightVars.find(sVariable) )
    {
      RooArgSet Row(*fSData->get(numEvent));
      totalYield += Row.getRealValue(sVariable); 
      
      return totalYield;
    }

  if( fSWeightVars.find(varname.c_str())  )
    {

      RooArgSet Row(*fSData->get(numEvent));
      totalYield += Row.getRealValue(varname.c_str() ); 
      
      return totalYield;
    }
  
  else
    coutE(InputArguments) << "InputVariable not in list of sWeighted variables" << endl;
    
  return -1;
}


//____________________________________________________________________
Double_t SPlot::GetSumOfEventSWeight(Int_t numEvent) const
{

  //Sum the SWeights for a particular event.
  //This sum should equal the total weight of that event.
  //This method is intended to be used as a check.


  if(numEvent > fSData->numEntries() )
    {
      coutE(InputArguments)  << "Invalid Entry Number" << endl;
      return -1;
    }

  if(numEvent < 0)
    {
      coutE(InputArguments)  << "Invalid Entry Number" << endl;
      return -1;    
    }

  Int_t numSWeightVars = this->GetNumSWeightVars();

  Double_t eventSWeight = 0;

  RooArgSet Row(*fSData->get(numEvent));
  
  for (Int_t i = 0; i < numSWeightVars; i++) 
    eventSWeight += Row.getRealValue(fSWeightVars.at(i)->GetName() );
 
  return  eventSWeight;
}


//_________________________________________________________________
Double_t SPlot::GetYieldFromSWeight(const char* sVariable) const
{


  //Sum the SWeights for a particular specie over all events
  //This should equal the total (weighted) yield of that specie
  //This method is intended as a check.

  Double_t totalYield = 0;

  std::string varname(sVariable);
  varname += "_sw";


  if(fSWeightVars.find(sVariable) )
    {
      for(Int_t i=0; i < fSData->numEntries(); i++)
	{
	  RooArgSet Row(*fSData->get(i));
	  totalYield += Row.getRealValue(sVariable); 
	}
      
      return totalYield;
    }

  if( fSWeightVars.find(varname.c_str())  )
    {
      for(Int_t i=0; i < fSData->numEntries(); i++)
	{
	  RooArgSet Row(*fSData->get(i));
	  totalYield += Row.getRealValue(varname.c_str() ); 
	}
      
      return totalYield;
    }
  
  else
    coutE(InputArguments) << "InputVariable not in list of sWeighted variables" << endl;
    
  return -1;
}


//______________________________________________________________________
RooArgList SPlot::GetSWeightVars() const
{
   
  //Return a RooArgList containing the SWeights

  RooArgList Args = fSWeightVars;
   
  return  Args;
   
}

//__________________________________________________________________ 
Int_t SPlot::GetNumSWeightVars() const
{
  //Return the number of SWeights
  //In other words, return the number of
  //species that we are trying to extract.
  
  RooArgList Args = fSWeightVars;
  
  return Args.getSize();
}


//____________________________________________________________________
void SPlot::AddSWeight( RooAbsPdf* pdf, const RooArgList &yieldsTmp, 
			const RooArgSet &projDeps, bool includeWeights) 
{  

  //
  // Method which adds the sWeights to the dataset.
  // Input is the PDF, a RooArgList of the yields (floating)
  // and a RooArgSet of the projDeps.
  //
  // The projDeps will not be normalized over when calculating the SWeights
  // and will be considered parameters, not observables.
  //
  // The SPlot will contain two new variables for each specie of name "varname":
  //
  // L_varname is the value of the pdf for the variable "varname" at values of this event
  // varname_sw is the value of the sWeight for the variable "varname" for this event
  
  // Find Parameters in the PDF to be considered fixed when calculating the SWeights
  // and be sure to NOT include the yields in that list
  //


  RooFit::MsgLevel currentLevel =  RooMsgService::instance().globalKillBelow();
  
  RooArgList* constParameters = (RooArgList*)pdf->getParameters(fSData) ;
  constParameters->remove(yieldsTmp, kTRUE, kTRUE);
  

  // Set these parameters constant and store them so they can later
  // be set to not constant
  std::vector<RooRealVar*> constVarHolder;  

  for(Int_t i = 0; i < constParameters->getSize(); i++) 
    {
      RooRealVar* varTemp = ( dynamic_cast<RooRealVar*>( constParameters->at(i) ) );
      if( varTemp->isConstant() == 0 )
	{
	  varTemp->setConstant();
	  constVarHolder.push_back(varTemp);
	}
    }
  
  // Fit yields to the data with all other variables held constant
  // This is necessary because SPlot assumes the yields minimixe -Log(likelihood)

  pdf->fitTo(*fSData, RooFit::Extended(kTRUE), RooFit::SumW2Error(kTRUE), RooFit::PrintLevel(-1), RooFit::PrintEvalErrors(-1) );

  // Hold the value of the fitted yields
  std::vector<double> yieldsHolder;
  
  for(Int_t i = 0; i < yieldsTmp.getSize(); i++)   
    yieldsHolder.push_back( ((RooRealVar*) yieldsTmp.at(i))->getVal());
  
  Int_t nspec = yieldsTmp.getSize();
  RooArgList yields = *(RooArgList*)yieldsTmp.snapshot(kFALSE);
  
  if(currentLevel <= RooFit::DEBUG)
    {
      coutI(InputArguments) << "Printing Yields" << endl;
      yields.Print();
    }
  
  // The list of variables to normalize over when calculating PDF values.

  RooArgSet vars(*fSData->get() );
  vars.remove(projDeps, kTRUE, kTRUE);
  
  // Attach data set

  // const_cast<RooAbsPdf*>(pdf)->attachDataSet(*fSData);

  pdf->attachDataSet(*fSData);

  // first calculate the pdf values for all species and all events
  std::vector<RooRealVar*> yieldvars ;
  RooArgSet* parameters = pdf->getParameters(fSData) ;

  std::vector<Double_t> yieldvalues ;
  for (Int_t k = 0; k < nspec; ++k) 
    {
      RooRealVar* thisyield = dynamic_cast<RooRealVar*>(yields.at(k)) ;
      RooRealVar* yieldinpdf = dynamic_cast<RooRealVar*>(parameters->find(thisyield->GetName() )) ;

      coutI(InputArguments)<< "yield in pdf: " << yieldinpdf->GetName() << " " << thisyield->getVal() << endl;
      
      yieldvars.push_back(yieldinpdf) ;
      yieldvalues.push_back(thisyield->getVal()) ;
    }
  
  Int_t numevents = fSData->numEntries() ;
  
  std::vector<std::vector<Double_t> > pdfvalues(numevents,std::vector<Double_t>(nspec,0)) ; 

      
  // set all yield to zero
  for(Int_t m=0; m<nspec; ++m) yieldvars[m]->setVal(0) ;
   

  // For every event and for every specie,
  // calculate the value of the component pdf for that specie
  // by setting the yield of that specie to 1
  // and all others to 0.  Evaluate the pdf for each event
  // and store the values.

  for (Int_t ievt = 0; ievt <numevents; ievt++) 
    {
      //   if (ievt % 100 == 0) 
      //  coutP(Eval)  << ".";


      //FIX THIS PART, EVALUATION PROGRESS!!

      RooStats::SetParameters(fSData->get(ievt), pdf->getVariables()); 

      //   RooArgSet row(*fSData->get(ievt));
       
      for(Int_t k = 0; k < nspec; ++k) 
	{
	  //Check that range of yields is at least (0,1), and fix otherwise
	  if(yieldvars[k]->getMin() > 0) 
	    {
	      coutW(InputArguments)  << "Minimum Range for " << yieldvars[k]->GetName() << " must be 0.  ";
	      coutW(InputArguments)  << "Setting min range to 0" << std::endl;
	      yieldvars[k]->setMin(0);
	    }

	  if(yieldvars[k]->getMax() < 1) 
	    {
	      coutW(InputArguments)  << "Maximum Range for " << yieldvars[k]->GetName() << " must be 1.  ";
	      coutW(InputArguments)  << "Setting max range to 1" << std::endl;
	      yieldvars[k]->setMax(1);
	    }

	  // set this yield to 1
	  yieldvars[k]->setVal( 1 ) ;
	  // evaluate the pdf    
	  Double_t f_k = pdf->getVal(&vars) ;
	  pdfvalues[ievt][k] = f_k ;
	  if( !(f_k>1 || f_k<1) ) 
	    coutW(InputArguments) << "Strange pdf value: " << ievt << " " << k << " " << f_k << std::endl ;
	  yieldvars[k]->setVal( 0 ) ;
	}
    }
   
  // check that the likelihood normalization is fine
  std::vector<Double_t> norm(nspec,0) ;
  for (Int_t ievt = 0; ievt <numevents ; ievt++) 
    {
      Double_t dnorm(0) ;
      for(Int_t k=0; k<nspec; ++k) dnorm += yieldvalues[k] * pdfvalues[ievt][k] ;
      for(Int_t j=0; j<nspec; ++j) norm[j] += pdfvalues[ievt][j]/dnorm ;
    }
   
  coutI(Contents) << "likelihood norms: "  ;

  for(Int_t k=0; k<nspec; ++k)  coutI(Contents) << norm[k] << " " ;
  coutI(Contents) << std::endl ;
   
  // Make a TMatrixD to hold the covariance matrix.
  TMatrixD covInv(nspec, nspec);
  for (Int_t i = 0; i < nspec; i++) for (Int_t j = 0; j < nspec; j++) covInv(i,j) = 0;
   
  coutI(Contents) << "Calculating covariance matrix";


  // Calculate the inverse covariance matrix, using weights
  for (Int_t ievt = 0; ievt < numevents; ++ievt) 
    {
      
      fSData->get(ievt) ;
     
      // Calculate contribution to the inverse of the covariance
      // matrix. See BAD 509 V2 eqn. 15

      // Sum for the denominator
      Double_t dsum(0);
      for(Int_t k = 0; k < nspec; ++k) 
	dsum += pdfvalues[ievt][k] * yieldvalues[k] ;
       
      for(Int_t n=0; n<nspec; ++n)
	for(Int_t j=0; j<nspec; ++j) 
	  {
	    if(includeWeights == kTRUE)
	      covInv(n,j) +=  fSData->weight()*pdfvalues[ievt][n]*pdfvalues[ievt][j]/(dsum*dsum) ;
	    else 
	      covInv(n,j) +=  pdfvalues[ievt][n]*pdfvalues[ievt][j]/(dsum*dsum) ;
	  }

      //ADDED WEIGHT ABOVE

    }
   
  // Covariance inverse should now be computed!
   
  // Invert to get the covariance matrix
  if (covInv.Determinant() <=0) 
    {
      coutE(Eval) << "SPlot Error: covariance matrix is singular; I can't invert it!" << std::endl;
      covInv.Print();
      return;
    }
   
  TMatrixD covMatrix(TMatrixD::kInverted,covInv);
   
  //check cov normalization
  if(currentLevel <= RooFit::DEBUG)
    {
      coutI(Eval) << "Checking Likelihood normalization:  " << std::endl;
      coutI(Eval) << "Yield of specie  Sum of Row in Matrix   Norm" << std::endl; 
      for(Int_t k=0; k<nspec; ++k) 
	{
	  Double_t covnorm(0) ;
	  for(Int_t m=0; m<nspec; ++m) covnorm += covInv[k][m]*yieldvalues[m] ;
	  Double_t sumrow(0) ;
	  for(Int_t m = 0; m < nspec; ++m) sumrow += covMatrix[k][m] ;
	  coutI(Eval)  << yieldvalues[k] << " " << sumrow << " " << covnorm << endl ;
	}
    }
  
  // calculate for each event the sWeight (BAD 509 V2 eq. 21)
  coutI(Eval) << "Calculating sWeight" << std::endl;
  std::vector<RooRealVar*> sweightvec ;
  std::vector<RooRealVar*> pdfvec ;  
  RooArgSet sweightset ;
   
  char wname[256] ;

  // Create and label the variables
  // used to store the SWeights

  fSWeightVars.Clear();

  for(Int_t k=0; k<nspec; ++k) 
    {
      sprintf(wname,"%s_sw", yieldvars[k]->GetName()) ;
      RooRealVar* var = new RooRealVar(wname,wname,0) ;
      sweightvec.push_back( var) ;
      sweightset.add(*var) ;
      fSWeightVars.add(*var);
    
      sprintf(wname,"L_%s", yieldvars[k]->GetName()) ;
      var = new RooRealVar(wname,wname,0) ;
      pdfvec.push_back( var) ;
      sweightset.add(*var) ;
    }

  // Create and fill a RooDataSet
  // with the SWeights
 
  RooDataSet* sWeightData = new RooDataSet("dataset", "dataset with sWeights", sweightset);
  
  for(Int_t ievt = 0; ievt < numevents; ++ievt) 
    {

      fSData->get(ievt) ;
       
      // sum for denominator
      Double_t dsum(0);
      for(Int_t k = 0; k < nspec; ++k)   dsum +=  pdfvalues[ievt][k] * yieldvalues[k] ;
      // covariance weighted pdf for each specief
      for(Int_t n=0; n<nspec; ++n) 
	{
	  Double_t nsum(0) ;
	  for(Int_t j=0; j<nspec; ++j) nsum += covMatrix(n,j) * pdfvalues[ievt][j] ;     
	   

	  //Add the sWeights here!!
	  //Include weights,
	  //ie events weights are absorbed into sWeight


	  if(includeWeights == kTRUE) sweightvec[n]->setVal(fSData->weight() * nsum/dsum) ;
	  else  sweightvec[n]->setVal( nsum/dsum) ;

	  pdfvec[n]->setVal( pdfvalues[ievt][n] ) ;

	  if( !(fabs(nsum/dsum)>=0 ) ) 
	    {
	      coutE(Contents) << "error: " << nsum/dsum << endl ;
	      return;
	    }
	}
      
      sWeightData->add(sweightset) ;
    }

    
  // Add the SWeights to the original data set

  fSData->merge(sWeightData);

  //Restore yield values

  for(Int_t i = 0; i < yieldsTmp.getSize(); i++) 
    ((RooRealVar*) yieldsTmp.at(i))->setVal(yieldsHolder.at(i));
   
  //Make any variables that were forced to constant no longer constant

  for(Int_t i=0; i < (Int_t) constVarHolder.size(); i++) 
    constVarHolder.at(i)->setConstant(kFALSE);

  return;

}
 SPlot.cxx:1
 SPlot.cxx:2
 SPlot.cxx:3
 SPlot.cxx:4
 SPlot.cxx:5
 SPlot.cxx:6
 SPlot.cxx:7
 SPlot.cxx:8
 SPlot.cxx:9
 SPlot.cxx:10
 SPlot.cxx:11
 SPlot.cxx:12
 SPlot.cxx:13
 SPlot.cxx:14
 SPlot.cxx:15
 SPlot.cxx:16
 SPlot.cxx:17
 SPlot.cxx:18
 SPlot.cxx:19
 SPlot.cxx:20
 SPlot.cxx:21
 SPlot.cxx:22
 SPlot.cxx:23
 SPlot.cxx:24
 SPlot.cxx:25
 SPlot.cxx:26
 SPlot.cxx:27
 SPlot.cxx:28
 SPlot.cxx:29
 SPlot.cxx:30
 SPlot.cxx:31
 SPlot.cxx:32
 SPlot.cxx:33
 SPlot.cxx:34
 SPlot.cxx:35
 SPlot.cxx:36
 SPlot.cxx:37
 SPlot.cxx:38
 SPlot.cxx:39
 SPlot.cxx:40
 SPlot.cxx:41
 SPlot.cxx:42
 SPlot.cxx:43
 SPlot.cxx:44
 SPlot.cxx:45
 SPlot.cxx:46
 SPlot.cxx:47
 SPlot.cxx:48
 SPlot.cxx:49
 SPlot.cxx:50
 SPlot.cxx:51
 SPlot.cxx:52
 SPlot.cxx:53
 SPlot.cxx:54
 SPlot.cxx:55
 SPlot.cxx:56
 SPlot.cxx:57
 SPlot.cxx:58
 SPlot.cxx:59
 SPlot.cxx:60
 SPlot.cxx:61
 SPlot.cxx:62
 SPlot.cxx:63
 SPlot.cxx:64
 SPlot.cxx:65
 SPlot.cxx:66
 SPlot.cxx:67
 SPlot.cxx:68
 SPlot.cxx:69
 SPlot.cxx:70
 SPlot.cxx:71
 SPlot.cxx:72
 SPlot.cxx:73
 SPlot.cxx:74
 SPlot.cxx:75
 SPlot.cxx:76
 SPlot.cxx:77
 SPlot.cxx:78
 SPlot.cxx:79
 SPlot.cxx:80
 SPlot.cxx:81
 SPlot.cxx:82
 SPlot.cxx:83
 SPlot.cxx:84
 SPlot.cxx:85
 SPlot.cxx:86
 SPlot.cxx:87
 SPlot.cxx:88
 SPlot.cxx:89
 SPlot.cxx:90
 SPlot.cxx:91
 SPlot.cxx:92
 SPlot.cxx:93
 SPlot.cxx:94
 SPlot.cxx:95
 SPlot.cxx:96
 SPlot.cxx:97
 SPlot.cxx:98
 SPlot.cxx:99
 SPlot.cxx:100
 SPlot.cxx:101
 SPlot.cxx:102
 SPlot.cxx:103
 SPlot.cxx:104
 SPlot.cxx:105
 SPlot.cxx:106
 SPlot.cxx:107
 SPlot.cxx:108
 SPlot.cxx:109
 SPlot.cxx:110
 SPlot.cxx:111
 SPlot.cxx:112
 SPlot.cxx:113
 SPlot.cxx:114
 SPlot.cxx:115
 SPlot.cxx:116
 SPlot.cxx:117
 SPlot.cxx:118
 SPlot.cxx:119
 SPlot.cxx:120
 SPlot.cxx:121
 SPlot.cxx:122
 SPlot.cxx:123
 SPlot.cxx:124
 SPlot.cxx:125
 SPlot.cxx:126
 SPlot.cxx:127
 SPlot.cxx:128
 SPlot.cxx:129
 SPlot.cxx:130
 SPlot.cxx:131
 SPlot.cxx:132
 SPlot.cxx:133
 SPlot.cxx:134
 SPlot.cxx:135
 SPlot.cxx:136
 SPlot.cxx:137
 SPlot.cxx:138
 SPlot.cxx:139
 SPlot.cxx:140
 SPlot.cxx:141
 SPlot.cxx:142
 SPlot.cxx:143
 SPlot.cxx:144
 SPlot.cxx:145
 SPlot.cxx:146
 SPlot.cxx:147
 SPlot.cxx:148
 SPlot.cxx:149
 SPlot.cxx:150
 SPlot.cxx:151
 SPlot.cxx:152
 SPlot.cxx:153
 SPlot.cxx:154
 SPlot.cxx:155
 SPlot.cxx:156
 SPlot.cxx:157
 SPlot.cxx:158
 SPlot.cxx:159
 SPlot.cxx:160
 SPlot.cxx:161
 SPlot.cxx:162
 SPlot.cxx:163
 SPlot.cxx:164
 SPlot.cxx:165
 SPlot.cxx:166
 SPlot.cxx:167
 SPlot.cxx:168
 SPlot.cxx:169
 SPlot.cxx:170
 SPlot.cxx:171
 SPlot.cxx:172
 SPlot.cxx:173
 SPlot.cxx:174
 SPlot.cxx:175
 SPlot.cxx:176
 SPlot.cxx:177
 SPlot.cxx:178
 SPlot.cxx:179
 SPlot.cxx:180
 SPlot.cxx:181
 SPlot.cxx:182
 SPlot.cxx:183
 SPlot.cxx:184
 SPlot.cxx:185
 SPlot.cxx:186
 SPlot.cxx:187
 SPlot.cxx:188
 SPlot.cxx:189
 SPlot.cxx:190
 SPlot.cxx:191
 SPlot.cxx:192
 SPlot.cxx:193
 SPlot.cxx:194
 SPlot.cxx:195
 SPlot.cxx:196
 SPlot.cxx:197
 SPlot.cxx:198
 SPlot.cxx:199
 SPlot.cxx:200
 SPlot.cxx:201
 SPlot.cxx:202
 SPlot.cxx:203
 SPlot.cxx:204
 SPlot.cxx:205
 SPlot.cxx:206
 SPlot.cxx:207
 SPlot.cxx:208
 SPlot.cxx:209
 SPlot.cxx:210
 SPlot.cxx:211
 SPlot.cxx:212
 SPlot.cxx:213
 SPlot.cxx:214
 SPlot.cxx:215
 SPlot.cxx:216
 SPlot.cxx:217
 SPlot.cxx:218
 SPlot.cxx:219
 SPlot.cxx:220
 SPlot.cxx:221
 SPlot.cxx:222
 SPlot.cxx:223
 SPlot.cxx:224
 SPlot.cxx:225
 SPlot.cxx:226
 SPlot.cxx:227
 SPlot.cxx:228
 SPlot.cxx:229
 SPlot.cxx:230
 SPlot.cxx:231
 SPlot.cxx:232
 SPlot.cxx:233
 SPlot.cxx:234
 SPlot.cxx:235
 SPlot.cxx:236
 SPlot.cxx:237
 SPlot.cxx:238
 SPlot.cxx:239
 SPlot.cxx:240
 SPlot.cxx:241
 SPlot.cxx:242
 SPlot.cxx:243
 SPlot.cxx:244
 SPlot.cxx:245
 SPlot.cxx:246
 SPlot.cxx:247
 SPlot.cxx:248
 SPlot.cxx:249
 SPlot.cxx:250
 SPlot.cxx:251
 SPlot.cxx:252
 SPlot.cxx:253
 SPlot.cxx:254
 SPlot.cxx:255
 SPlot.cxx:256
 SPlot.cxx:257
 SPlot.cxx:258
 SPlot.cxx:259
 SPlot.cxx:260
 SPlot.cxx:261
 SPlot.cxx:262
 SPlot.cxx:263
 SPlot.cxx:264
 SPlot.cxx:265
 SPlot.cxx:266
 SPlot.cxx:267
 SPlot.cxx:268
 SPlot.cxx:269
 SPlot.cxx:270
 SPlot.cxx:271
 SPlot.cxx:272
 SPlot.cxx:273
 SPlot.cxx:274
 SPlot.cxx:275
 SPlot.cxx:276
 SPlot.cxx:277
 SPlot.cxx:278
 SPlot.cxx:279
 SPlot.cxx:280
 SPlot.cxx:281
 SPlot.cxx:282
 SPlot.cxx:283
 SPlot.cxx:284
 SPlot.cxx:285
 SPlot.cxx:286
 SPlot.cxx:287
 SPlot.cxx:288
 SPlot.cxx:289
 SPlot.cxx:290
 SPlot.cxx:291
 SPlot.cxx:292
 SPlot.cxx:293
 SPlot.cxx:294
 SPlot.cxx:295
 SPlot.cxx:296
 SPlot.cxx:297
 SPlot.cxx:298
 SPlot.cxx:299
 SPlot.cxx:300
 SPlot.cxx:301
 SPlot.cxx:302
 SPlot.cxx:303
 SPlot.cxx:304
 SPlot.cxx:305
 SPlot.cxx:306
 SPlot.cxx:307
 SPlot.cxx:308
 SPlot.cxx:309
 SPlot.cxx:310
 SPlot.cxx:311
 SPlot.cxx:312
 SPlot.cxx:313
 SPlot.cxx:314
 SPlot.cxx:315
 SPlot.cxx:316
 SPlot.cxx:317
 SPlot.cxx:318
 SPlot.cxx:319
 SPlot.cxx:320
 SPlot.cxx:321
 SPlot.cxx:322
 SPlot.cxx:323
 SPlot.cxx:324
 SPlot.cxx:325
 SPlot.cxx:326
 SPlot.cxx:327
 SPlot.cxx:328
 SPlot.cxx:329
 SPlot.cxx:330
 SPlot.cxx:331
 SPlot.cxx:332
 SPlot.cxx:333
 SPlot.cxx:334
 SPlot.cxx:335
 SPlot.cxx:336
 SPlot.cxx:337
 SPlot.cxx:338
 SPlot.cxx:339
 SPlot.cxx:340
 SPlot.cxx:341
 SPlot.cxx:342
 SPlot.cxx:343
 SPlot.cxx:344
 SPlot.cxx:345
 SPlot.cxx:346
 SPlot.cxx:347
 SPlot.cxx:348
 SPlot.cxx:349
 SPlot.cxx:350
 SPlot.cxx:351
 SPlot.cxx:352
 SPlot.cxx:353
 SPlot.cxx:354
 SPlot.cxx:355
 SPlot.cxx:356
 SPlot.cxx:357
 SPlot.cxx:358
 SPlot.cxx:359
 SPlot.cxx:360
 SPlot.cxx:361
 SPlot.cxx:362
 SPlot.cxx:363
 SPlot.cxx:364
 SPlot.cxx:365
 SPlot.cxx:366
 SPlot.cxx:367
 SPlot.cxx:368
 SPlot.cxx:369
 SPlot.cxx:370
 SPlot.cxx:371
 SPlot.cxx:372
 SPlot.cxx:373
 SPlot.cxx:374
 SPlot.cxx:375
 SPlot.cxx:376
 SPlot.cxx:377
 SPlot.cxx:378
 SPlot.cxx:379
 SPlot.cxx:380
 SPlot.cxx:381
 SPlot.cxx:382
 SPlot.cxx:383
 SPlot.cxx:384
 SPlot.cxx:385
 SPlot.cxx:386
 SPlot.cxx:387
 SPlot.cxx:388
 SPlot.cxx:389
 SPlot.cxx:390
 SPlot.cxx:391
 SPlot.cxx:392
 SPlot.cxx:393
 SPlot.cxx:394
 SPlot.cxx:395
 SPlot.cxx:396
 SPlot.cxx:397
 SPlot.cxx:398
 SPlot.cxx:399
 SPlot.cxx:400
 SPlot.cxx:401
 SPlot.cxx:402
 SPlot.cxx:403
 SPlot.cxx:404
 SPlot.cxx:405
 SPlot.cxx:406
 SPlot.cxx:407
 SPlot.cxx:408
 SPlot.cxx:409
 SPlot.cxx:410
 SPlot.cxx:411
 SPlot.cxx:412
 SPlot.cxx:413
 SPlot.cxx:414
 SPlot.cxx:415
 SPlot.cxx:416
 SPlot.cxx:417
 SPlot.cxx:418
 SPlot.cxx:419
 SPlot.cxx:420
 SPlot.cxx:421
 SPlot.cxx:422
 SPlot.cxx:423
 SPlot.cxx:424
 SPlot.cxx:425
 SPlot.cxx:426
 SPlot.cxx:427
 SPlot.cxx:428
 SPlot.cxx:429
 SPlot.cxx:430
 SPlot.cxx:431
 SPlot.cxx:432
 SPlot.cxx:433
 SPlot.cxx:434
 SPlot.cxx:435
 SPlot.cxx:436
 SPlot.cxx:437
 SPlot.cxx:438
 SPlot.cxx:439
 SPlot.cxx:440
 SPlot.cxx:441
 SPlot.cxx:442
 SPlot.cxx:443
 SPlot.cxx:444
 SPlot.cxx:445
 SPlot.cxx:446
 SPlot.cxx:447
 SPlot.cxx:448
 SPlot.cxx:449
 SPlot.cxx:450
 SPlot.cxx:451
 SPlot.cxx:452
 SPlot.cxx:453
 SPlot.cxx:454
 SPlot.cxx:455
 SPlot.cxx:456
 SPlot.cxx:457
 SPlot.cxx:458
 SPlot.cxx:459
 SPlot.cxx:460
 SPlot.cxx:461
 SPlot.cxx:462
 SPlot.cxx:463
 SPlot.cxx:464
 SPlot.cxx:465
 SPlot.cxx:466
 SPlot.cxx:467
 SPlot.cxx:468
 SPlot.cxx:469
 SPlot.cxx:470
 SPlot.cxx:471
 SPlot.cxx:472
 SPlot.cxx:473
 SPlot.cxx:474
 SPlot.cxx:475
 SPlot.cxx:476
 SPlot.cxx:477
 SPlot.cxx:478
 SPlot.cxx:479
 SPlot.cxx:480
 SPlot.cxx:481
 SPlot.cxx:482
 SPlot.cxx:483
 SPlot.cxx:484
 SPlot.cxx:485
 SPlot.cxx:486
 SPlot.cxx:487
 SPlot.cxx:488
 SPlot.cxx:489
 SPlot.cxx:490
 SPlot.cxx:491
 SPlot.cxx:492
 SPlot.cxx:493
 SPlot.cxx:494
 SPlot.cxx:495
 SPlot.cxx:496
 SPlot.cxx:497
 SPlot.cxx:498
 SPlot.cxx:499
 SPlot.cxx:500
 SPlot.cxx:501
 SPlot.cxx:502
 SPlot.cxx:503
 SPlot.cxx:504
 SPlot.cxx:505
 SPlot.cxx:506
 SPlot.cxx:507
 SPlot.cxx:508
 SPlot.cxx:509
 SPlot.cxx:510
 SPlot.cxx:511
 SPlot.cxx:512
 SPlot.cxx:513
 SPlot.cxx:514
 SPlot.cxx:515
 SPlot.cxx:516
 SPlot.cxx:517
 SPlot.cxx:518
 SPlot.cxx:519
 SPlot.cxx:520
 SPlot.cxx:521
 SPlot.cxx:522
 SPlot.cxx:523
 SPlot.cxx:524
 SPlot.cxx:525
 SPlot.cxx:526
 SPlot.cxx:527
 SPlot.cxx:528
 SPlot.cxx:529
 SPlot.cxx:530
 SPlot.cxx:531
 SPlot.cxx:532
 SPlot.cxx:533
 SPlot.cxx:534
 SPlot.cxx:535
 SPlot.cxx:536
 SPlot.cxx:537
 SPlot.cxx:538
 SPlot.cxx:539
 SPlot.cxx:540
 SPlot.cxx:541
 SPlot.cxx:542
 SPlot.cxx:543
 SPlot.cxx:544
 SPlot.cxx:545
 SPlot.cxx:546
 SPlot.cxx:547
 SPlot.cxx:548
 SPlot.cxx:549
 SPlot.cxx:550
 SPlot.cxx:551
 SPlot.cxx:552
 SPlot.cxx:553
 SPlot.cxx:554
 SPlot.cxx:555
 SPlot.cxx:556
 SPlot.cxx:557
 SPlot.cxx:558
 SPlot.cxx:559
 SPlot.cxx:560
 SPlot.cxx:561
 SPlot.cxx:562
 SPlot.cxx:563
 SPlot.cxx:564
 SPlot.cxx:565
 SPlot.cxx:566
 SPlot.cxx:567
 SPlot.cxx:568
 SPlot.cxx:569
 SPlot.cxx:570
 SPlot.cxx:571
 SPlot.cxx:572
 SPlot.cxx:573
 SPlot.cxx:574
 SPlot.cxx:575
 SPlot.cxx:576
 SPlot.cxx:577
 SPlot.cxx:578
 SPlot.cxx:579
 SPlot.cxx:580
 SPlot.cxx:581
 SPlot.cxx:582
 SPlot.cxx:583
 SPlot.cxx:584
 SPlot.cxx:585
 SPlot.cxx:586
 SPlot.cxx:587
 SPlot.cxx:588
 SPlot.cxx:589
 SPlot.cxx:590
 SPlot.cxx:591
 SPlot.cxx:592
 SPlot.cxx:593
 SPlot.cxx:594
 SPlot.cxx:595
 SPlot.cxx:596
 SPlot.cxx:597
 SPlot.cxx:598
 SPlot.cxx:599
 SPlot.cxx:600
 SPlot.cxx:601
 SPlot.cxx:602
 SPlot.cxx:603
 SPlot.cxx:604
 SPlot.cxx:605
 SPlot.cxx:606
 SPlot.cxx:607
 SPlot.cxx:608
 SPlot.cxx:609
 SPlot.cxx:610
 SPlot.cxx:611
 SPlot.cxx:612
 SPlot.cxx:613
 SPlot.cxx:614
 SPlot.cxx:615
 SPlot.cxx:616
 SPlot.cxx:617
 SPlot.cxx:618
 SPlot.cxx:619
 SPlot.cxx:620
 SPlot.cxx:621
 SPlot.cxx:622
 SPlot.cxx:623
 SPlot.cxx:624