C++ DataPoint::GetObservable方法代码示例

//Return a list of data points
//Each should take the data average value of each continuous observable
//Each should represent one combination of possible discrete values
vector<DataPoint*> PhaseSpaceBoundary::GetDiscreteCombinations() const
{
	if( storedCombinationID == uniqueID && !StoredCombinations.empty() )
	{
		return StoredCombinations;
	}

	while( !StoredCombinations.empty() )
	{
		if( StoredCombinations.back() != NULL ) delete StoredCombinations.back();
		StoredCombinations.pop_back();
	}

	//Calculate all possible combinations of discrete observables
	vector<string> thisAllNames = this->GetAllNames();
	vector<vector<double> > discreteValues;
	vector<string> discreteNames, continuousNames;
	vector<vector<double> > discreteCombinations = StatisticsFunctions::DiscreteCombinations( &thisAllNames, this, discreteNames, continuousNames, discreteValues );

	(void) continuousNames;

	//Create the data points to return
	vector<DataPoint*> newDataPoints;

	DataPoint* tempPoint = this->GetMidPoint();
	if( tempPoint == NULL ) return newDataPoints;

	for( unsigned int combinationIndex = 0; combinationIndex < discreteCombinations.size(); ++combinationIndex )
	{
		DataPoint* templateDataPoint = new DataPoint( *tempPoint );

		//Output the discrete values for this combination
		for( unsigned int discreteIndex = 0; discreteIndex < discreteNames.size(); ++discreteIndex )
		{
			//Set the data point
			Observable* oldValue = templateDataPoint->GetObservable( discreteNames[discreteIndex] );

			Observable* newValue = new Observable( oldValue->GetName(), discreteCombinations[combinationIndex][discreteIndex], oldValue->GetUnit() );

			templateDataPoint->SetObservable( discreteNames[discreteIndex], newValue );
			delete newValue;
		}

		newDataPoints.push_back( templateDataPoint );
	}

	delete tempPoint;

	StoredCombinations = newDataPoints;

	storedCombinationID = uniqueID; 

	return newDataPoints;
}

vector<DataPoint*> RapidFitIntegrator::getGSLIntegrationPoints( unsigned int number, vector<double> maxima, vector<double> minima, DataPoint* templateDataPoint, vector<string> doIntegrate, PhaseSpaceBoundary* thisBound )
{
	pthread_mutex_lock( &GSL_DATAPOINT_GET_THREADLOCK );

	bool pointsGood = true;
	if( !_global_doEval_points.empty() )
	{
		if( _global_doEval_points[0] != NULL )
		{
			vector<string> allConsts = thisBound->GetDiscreteNames();
			for( unsigned int i=0; i< allConsts.size(); ++i )
			{
				double haveVal = _global_doEval_points[0]->GetObservable( allConsts[i] )->GetValue();
				double wantVal = templateDataPoint->GetObservable( allConsts[i] )->GetValue();

				if( haveVal != wantVal )
				{
					pointsGood = false;
					break;
				}
			}
		}
	}

	if( ( number != _global_doEval_points.size() ) || (( ( _global_range_minima != minima ) || ( _global_range_maxima != maxima ) ) || !pointsGood ) )
	{
		clearGSLIntegrationPoints();
		_global_doEval_points = initGSLDataPoints( number, maxima, minima, templateDataPoint, doIntegrate );
		_global_range_minima = minima;
		_global_range_maxima = maxima;
		_global_observable_names = doIntegrate;
	}

	vector<string> allObs = _global_doEval_points[0]->GetAllNames();
	for( unsigned int i=0; i< _global_doEval_points.size(); ++i )
	{
		DataPoint* thisPoint = _global_doEval_points[i];
		for( unsigned int j=0; j< allObs.size(); ++j )
		{
			Observable* thisObs = thisPoint->GetObservable( j );
			thisObs->SetBinNumber( -1 );
			thisObs->SetBkgBinNumber( -1 );
		}
		thisPoint->ClearPerEventData();
	}

	pthread_mutex_unlock( &GSL_DATAPOINT_GET_THREADLOCK );
	return _global_doEval_points;
}

void LegendreMomentShape::Generate(IDataSet* dataSet, const PhaseSpaceBoundary* boundary, const std::string mKKname, const std::string phiname, const std::string ctheta_1name, const std::string ctheta_2name)
{
	double**** c    = newcoefficients();
	double**** c_sq = newcoefficients(); // Used in caclulating the error
	mKK_min = boundary->GetConstraint(mKKname)->GetMinimum();
	mKK_max = boundary->GetConstraint(mKKname)->GetMaximum();
	const double mK  = 0.493677; // TODO: read these from config somehow
	const double mBs  = 5.36677;
	const double mPhi= 1.019461;
	int numEvents = dataSet->GetDataNumber();
	// Calculate the coefficients by summing over the dataset
	std::cout << "Sum over " << numEvents << " events" << std::endl;
	for (int e = 0; e < numEvents; e++)
	{
		// Retrieve the data point
		DataPoint * event = dataSet->GetDataPoint(e);
		double phi        = event->GetObservable(phiname)->GetValue();
		double ctheta_1   = event->GetObservable(ctheta_1name)->GetValue();
		double ctheta_2   = event->GetObservable(ctheta_2name)->GetValue();
		double mKK        = event->GetObservable(mKKname)->GetValue();
		double mKK_mapped = (mKK - mKK_min)/(mKK_max-mKK_min)*2.+ (-1);
		// Calculate phase space element
		double p1_st = DPHelpers::daughterMomentum(mKK, mK, mK);
		double p3    = DPHelpers::daughterMomentum(mBs,mKK,mPhi);
		double val = p1_st*p3;
		for ( int l = 0; l < l_max; l++ )
			for ( int i = 0; i < i_max; i++ )
				for ( int k = 0; k < k_max; k++ )
					for ( int j = 0; j < j_max; j++ )
					{
						double coeff = ((2*l + 1)/2.)*((2*i + 1)/2.)*Moment(l, i, k, j, mKK_mapped, phi, ctheta_1, ctheta_2)/val;
						c[l][i][k][j] += coeff;
						c_sq[l][i][k][j] += coeff * coeff;
					}
	}
	// Accept or reject the coefficients
	double threshold = 4; // TODO: read from config
	std::cout << "Keeping coefficients more significant than " << threshold << "σ" << std::endl;
	for ( int l = 0; l < l_max; l++ )
		for ( int i = 0; i < i_max; i++ )
			for ( int k = 0; k < k_max; k++ )
				for ( int j = 0; j < j_max; j++ )
				{
					if(std::abs(c[l][i][k][j]) < 1e-12)
					{
						c[l][i][k][j] = 0.;
						continue;
					}
					double error = sqrt(1./numEvents/numEvents * ( c_sq[l][i][k][j] - c[l][i][k][j]*c[l][i][k][j]/numEvents) );
					double signif = std::abs(c[l][i][k][j]/numEvents)/error;
					if ( signif > threshold )
					{
						printf("c[%d][%d][%d][%d] = %f;// ± %f with significance %fσ\n", l, i, k, j, c[l][i][k][j]/numEvents, error, signif );
						c[l][i][k][j] /= numEvents;
					}
					else
						c[l][i][k][j] = 0.;
				}
	storecoefficients(c);
	deletecoefficients(c);
	deletecoefficients(c_sq);
	init = false;
}

vector<DataPoint*> RapidFitIntegrator::initGSLDataPoints( unsigned int number, vector<double> maxima, vector<double> minima, DataPoint* templateDataPoint, vector<string> doIntegrate )
{
	vector<DataPoint*> doEval_points;
	pthread_mutex_lock( &GSL_DATAPOINT_MANAGEMENT_THREADLOCK );

#ifdef __RAPIDFIT_USE_GSL

	unsigned int nDim = (unsigned) minima.size();

	unsigned int npoint = number;
	vector<double> * integrationPoints = new vector<double>[ nDim ];

	//pthread_mutex_lock( &gsl_mutex );
	//cout << "Allocating GSL Integration Tool. nDim " << nDim << endl;
	gsl_qrng * q = NULL;
	try
	{
		//q = gsl_qrng_alloc( gsl_qrng_sobol, (unsigned)nDim );
		q = gsl_qrng_alloc( gsl_qrng_niederreiter_2, (unsigned)nDim );
	}
	catch(...)
	{
		cout << "Can't Allocate Integration Tool for GSL Integral." << endl;
		cout << " Dim: " << nDim << endl;
		exit(-742);
	}

	if( q == NULL )
	{
		cout << "Can't Allocate Integration Tool for GSL Integral." << endl;
		cout << " Dim: " << nDim << endl;
		exit(-741);
	}

	double* v = new double[ nDim ];
	for( unsigned int i = 0; i < npoint; i++)
	{
		gsl_qrng_get( q, v );
		for( unsigned int j = 0; j < nDim; j++)
		{
			integrationPoints[j].push_back( v[j] );
		}
	}
	delete v;
	gsl_qrng_free(q);
	//cout << "Freed GSL Integration Tool" << endl;
	//pthread_mutex_unlock( &gsl_mutex );

	/*vector<double> minima_v, maxima_v;
	  for( unsigned int i=0; i< nDim; ++i )
	  {
	  minima_v.push_back( minima[i] );
	  maxima_v.push_back( maxima[i] );
	  }*/
	//cout << "Constructing Functions" << endl;
	//IntegratorFunction* quickFunction = new IntegratorFunction( functionToWrap, NewDataPoint, doIntegrate, dontIntegrate, NewBoundary, componentIndex, minima_v, maxima_v );

	for (unsigned int i = 0; i < integrationPoints[0].size(); ++i)
	{
		double* point = new double[ nDim ];
		for ( unsigned int j = 0; j < nDim; j++)
		{
			//cout << doIntegrate[j] << " " << maxima[j] << " " << minima[j] << " " << integrationPoints[j][i] << endl;
			point[j] = integrationPoints[j][i]*(maxima[j]-minima[j])+minima[j];
		}

		DataPoint* thisPoint = new DataPoint( *templateDataPoint );

		thisPoint->ClearPerEventData();
		vector<string> allObs=thisPoint->GetAllNames();
		for( unsigned int j=0; j<allObs.size(); ++j )
		{
			Observable* thisObs = thisPoint->GetObservable( j );
			thisObs->SetBinNumber(-1);
			thisObs->SetBkgBinNumber(-1);
		}

		for( unsigned int k=0; k< doIntegrate.size(); ++k )
		{
			thisPoint->SetObservable( doIntegrate[k], point[k], "noUnitsHere" );
		}

		doEval_points.push_back( thisPoint );
		//result += quickFunction->DoEval( point );
		delete[] point;
	}

	delete[] integrationPoints;
#endif

	(void) maxima; (void) minima; (void) number;
	pthread_mutex_unlock( &GSL_DATAPOINT_MANAGEMENT_THREADLOCK );
	return doEval_points;
}

//Constructor with correct argument
MakeFoam::MakeFoam( IPDF * InputPDF, PhaseSpaceBoundary * InputBoundary, DataPoint * InputPoint ) : finishedCells(), centerPoints(), centerValues(), cellIntegrals(), integratePDF(InputPDF)
{
	//Make the container to hold the possible cells
	queue<PhaseSpaceBoundary*> possibleCells;
	PhaseSpaceBoundary* firstCell = InputBoundary;
	possibleCells.push(firstCell);

	//Make a list of observables to integrate over
	vector<string> doIntegrate, dontIntegrate;
	vector<string> pdfDontIntegrate = InputPDF->GetDoNotIntegrateList();
	StatisticsFunctions::DoDontIntegrateLists( InputPDF, InputBoundary, &(pdfDontIntegrate), doIntegrate, dontIntegrate );

	//Continue until all possible cells have been examined
	while ( !possibleCells.empty() )
	{
		//Remove the next possible cell
		PhaseSpaceBoundary* currentCell = possibleCells.front();
		possibleCells.pop();

		//Set up the histogram storage
		vector< vector<double> > histogramBinHeights, histogramBinMiddles, histogramBinMaxes;
		double normalisation = 0.0;
		for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
		{
			vector<double> binHeights, binMiddles, binMaxes;
			IConstraint * temporaryConstraint = currentCell->GetConstraint( doIntegrate[observableIndex] );
			double minimum = temporaryConstraint->GetMinimum();
			double delta = ( temporaryConstraint->GetMaximum() - minimum ) / (double)HISTOGRAM_BINS;

			//Loop over bins
			for (int binIndex = 0; binIndex < HISTOGRAM_BINS; ++binIndex )
			{
				binHeights.push_back(0.0);
				binMiddles.push_back( minimum + ( delta * ( binIndex + 0.5 ) ) );
				binMaxes.push_back( minimum + ( delta * ( binIndex + 1.0 ) ) );
			}

			histogramBinHeights.push_back(binHeights);
			histogramBinMiddles.push_back(binMiddles);
			histogramBinMaxes.push_back(binMaxes);
		}

		//MC sample the cell, make projections, sort of
		for (int sampleIndex = 0; sampleIndex < MAXIMUM_SAMPLES; ++sampleIndex )
		{
			//Create a data point within the current cell
			DataPoint samplePoint( InputPoint->GetAllNames() );
			for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
			{
				//Generate random values to explore integrable observables
				IConstraint * temporaryConstraint = currentCell->GetConstraint( doIntegrate[observableIndex] );
				samplePoint.SetObservable( doIntegrate[observableIndex], temporaryConstraint->CreateObservable() );
			}
			for (unsigned int observableIndex = 0; observableIndex < dontIntegrate.size(); ++observableIndex )
			{
				//Use given values for unintegrable observables
				Observable * temporaryObservable = new Observable( *(InputPoint->GetObservable( dontIntegrate[observableIndex] )) );
				samplePoint.SetObservable( dontIntegrate[observableIndex], temporaryObservable );
				delete temporaryObservable;
			}

			//Evaluate the function at this point
			double sampleValue = InputPDF->Evaluate( &samplePoint );
			normalisation += sampleValue;

			//Populate the histogram
			for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
			{
				double observableValue = samplePoint.GetObservable( doIntegrate[observableIndex] )->GetValue();

				for ( int binIndex = 0; binIndex < HISTOGRAM_BINS; ++binIndex )
				{
					if ( observableValue < histogramBinMaxes[observableIndex][unsigned(binIndex)] )
					{
						histogramBinHeights[observableIndex][unsigned(binIndex)] += sampleValue;
						break;
					}
				}
			}
		}

		//Normalise the histograms
		for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
		{
			for ( int binIndex = 0; binIndex < HISTOGRAM_BINS; ++binIndex )
			{
				histogramBinHeights[observableIndex][unsigned(binIndex)] /= normalisation;
			}
		}

		//Find the maximum gradient
		vector<double> midPoints;
		string maximumGradientObservable, unit;
		double maximumGradient=0.;
		double lowPoint=0.;
		double splitPoint=0.;
		double highPoint=0.;
		for (unsigned int observableIndex = 0; observableIndex < doIntegrate.size(); ++observableIndex )
		{
//.........这里部分代码省略.........