c_cpp 用于阅读我和其他几个人使用的光谱格式的快速代码

Posted 2021-05-14

#ifndef DRAW_SINGLE_SPECTRA_C
#define DRAW_SINGLE_SPECTRA_C


#include "SpectraLoader.h"

#include <sys/stat.h>

 bool file_exists (const std::string& name) {
  struct stat buffer;   
  return (stat (name.c_str(), &buffer) == 0); 
}

string file_name( string energy, string plc, string charge, string iCen ){
	string base = "/Users/danielbrandenburg/bnl/local/data/RcpAnalysis/spectra/";
	return base + "spectra_" + energy + "_" + plc + "_" + charge + "_" + iCen + ".dat";

}

int nPtBins = 26;
double ptBins[] = { 
0.0,
0.5,
0.6, 
0.7, 
0.8, 
0.9, 
1.0, 
1.1, 
1.2, 
1.3, 
1.4, 
1.5, 
1.6, 
1.7, 
1.8, 
1.9, 
2.0	, 
2.2	, 
2.4	, 
2.6	, 
2.8	,
3.0 , 
3.5, 
4.5, 
5.0, 
6.0, 
6.8 };

TH1 * normalize_binning( TH1 * in ){
	
	string name = in->GetName();
	name += + "_normed";
	TH1 * out = new TH1D( name.c_str(), "", nPtBins, ptBins );


	DEBUG( "Input has " << in->GetNbinsX() << " bins" );
	for ( int i = 1; i <= in->GetNbinsX(); i++ ){
		out->SetBinContent( i + 1, in->GetBinContent( i ) );
		out->SetBinError( i + 1, in->GetBinError( i ) );
	}

	return out;
}


TH1* draw_single_spectra( 	string energy, string plc, string charge, string iCen,
							int color = kRed, string draw_opt = "", double scaler = 1.0 ){
	
	gStyle->SetOptStat( 0 );


	string fn = file_name( energy, plc, charge, iCen );
	if ( !file_exists( fn ) )
		return new TH1D( "err", "", nPtBins, ptBins );

	SpectraLoader sl( fn );

	TH1* sys = normalize_binning( sl.sysHisto( fn + "_sys" ) );
	TH1* stat = normalize_binning( sl.statHisto( fn + "_stat" ) );

	sys->Scale( scaler );
	stat->Scale( scaler );

	
	sys->SetTitle( " ; pT [GeV/c]; dN^{2} / ( N_{evt} 2 #pi pT dpT dy )" );
	
	sys->SetLineColor( color );
	sys->SetMarkerStyle( 8 );
	sys->SetMarkerColor( color );
	sys->SetFillColorAlpha( color, 0.5 );

	stat->SetLineColor( color );
	stat->SetMarkerStyle( 8 );
	stat->SetMarkerColor( color );

	sys->Draw( (draw_opt + " e2").c_str() );
	stat->Draw( "same" );

	gPad->SetLogy(1);

	return sys;
}

#endif

#include <fstream>
#include "TGraphErrors.h"

class SpectraLoader{

public:
	vector<double> pT;
	vector<double> value;
	vector<double> stat;
	vector<double> sys;

	// computed
	vector<double> width;

	SpectraLoader( string fname ){

		ifstream inf( fname.c_str() );

		string tmp;
		double a, b, c, d;

		getline( inf, tmp );

		while( inf >> a >> b >> c >> d ){

			pT.push_back( a );
			value.push_back( b );
			stat.push_back( c );
			sys.push_back( d );

		}	

		if ( pT[ 0 ] > pT[ pT.size() - 1 ] ){
			reverse( pT.begin(), pT.end() );
			reverse( value.begin(), value.end() );
			reverse( stat.begin(), stat.end() );
			reverse( sys.begin(), sys.end() );
		}


		for ( int i = 0; i < pT.size(); i++ ){

			if ( i < pT.size() - 1  )
				width.push_back( (pT[ i+1 ] - pT[ i ]) / 2.0 ) ;
			else
				width.push_back( (pT[ i ] - pT[ i - 1 ]) / 2.0 );
		}


		inf.close();
	}

	void trim( int N = 1 ){

		for ( int i = 0; i < N; i++ ){
			pT.erase( pT.begin() );
			width.erase( width.begin() );
			value.erase( value.begin() );
			stat.erase( stat.begin() );
			sys.erase( sys.begin() );
		}
	}

	void trunc( int N = 1 ){

		for ( int i = 0; i < N; i++ ){
			pT.erase( pT.end() );
			width.erase( width.end() );
			value.erase( value.end() );
			stat.erase( stat.end() );
			sys.erase( sys.end() );
		}
	}

	TGraphErrors * statGraph(){
		TGraphErrors * graph = new TGraphErrors( pT.size(), pT.data(), value.data(), width.data(), stat.data() );
		return graph;
	}

	vector<double> getBins() {
		vector<double> bins;
		// make the bin edges
		for ( int i = 0; i < pT.size(); i++ ){
			if ( 0 == i){
				bins.push_back( pT[ 0 ] - width[ 0 ] );
			}
			bins.push_back( pT[ i ] + width[ i ] );
		}

		for ( int i = 0; i < bins.size(); i++ ){
			INFO( tag, "bin edge [" << i << "] = " << bins[ i ] );
		}
		return bins;
	}

	TH1D * statHisto( string name ){
		vector<double> bins = getBins();
		TH1D * h = new TH1D( name.c_str(), "", bins.size() - 1, bins.data() );
		for ( int i = 0; i < bins.size() + 1; i++ ){
			h->SetBinContent( i + 1, value[ i ] );
			h->SetBinError( i + 1, stat[ i ] );
		}
		return h;
	}

	

	TH1D * sysHisto( string name ){
		vector<double> bins = getBins();
		TH1D * h = new TH1D( name.c_str(), "", bins.size() - 1, bins.data() );
		for ( int i = 0; i < bins.size() + 1; i++ ){
			h->SetBinContent( i + 1, value[ i ] );
			h->SetBinError( i + 1, sys[ i ] );
		}
		return h;
	}

	TH1D * relSysHisto( string name ){
		vector<double> bins = getBins();
		TH1D * h = new TH1D( name.c_str(), "", bins.size() - 1, bins.data() );
		for ( int i = 0; i < bins.size() + 1; i++ ){
			h->SetBinContent( i + 1, sys[i] / value[ i ] );
			h->SetBinError( i + 1, 0 );
		}
		return h;
	}

	TH1D * relStatHisto( string name ){
		vector<double> bins = getBins();
		TH1D * h = new TH1D( name.c_str(), "", bins.size() - 1, bins.data() );
		for ( int i = 0; i < bins.size() + 1; i++ ){
			h->SetBinContent( i + 1, stat[i] / value[ i ] );
			h->SetBinError( i + 1, 0 );
		}
		return h;
	}

};