File indexing completed on 2024-09-27 07:03:41
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036 #include <iostream>
0037 #include <fstream>
0038
0039 #include "readinluminosity.h"
0040 #include "starlightconstants.h"
0041 #include "inputParameters.h"
0042
0043 using namespace std;
0044
0045
0046
0047 readLuminosity::readLuminosity(const inputParameters& inputParametersInstance)
0048 : _Warray(0), _Yarray(0), _Farray(0), _Farray1(0), _Farray2(0),
0049 _f_WYarray(0), _g_Earray(0)
0050 , _ReadInputNPT(inputParametersInstance.nmbPtBinsInterference())
0051 , _ReadInputnumy(inputParametersInstance.nmbRapidityBins())
0052 , _ReadInputnumw(inputParametersInstance.nmbWBins())
0053 , _ReadInputnumega(inputParametersInstance.nmbEnergyBins())
0054 , _ReadInputnumQ2(inputParametersInstance.nmbGammaQ2Bins())
0055 , _ReadInputgg_or_gP(inputParametersInstance.productionMode())
0056 , _ReadInputinterferencemode(inputParametersInstance.interferenceEnabled())
0057 , _baseFileName(inputParametersInstance.baseFileName())
0058 {
0059
0060 }
0061
0062
0063
0064 readLuminosity::~readLuminosity()
0065 {
0066 if(_Warray) delete [] _Warray;
0067 if(_Yarray) delete [] _Yarray;
0068 if(_Farray) delete [] _Farray;
0069 if(_Farray1) delete [] _Farray1;
0070 if(_Farray2) delete [] _Farray2;
0071
0072 if(_f_WYarray) delete [] _f_WYarray;
0073 if(_g_Earray) delete [] _g_Earray;
0074 if(_g_EQ2array) delete _g_EQ2array;
0075 }
0076
0077
0078
0079 void readLuminosity::read()
0080 {
0081
0082 if(!_Warray) _Warray = new double[_ReadInputnumw];
0083 if(!_Yarray) _Yarray = new double[_ReadInputnumy];
0084 if(!_Farray)
0085 {
0086 _Farray = new double*[_ReadInputnumw];
0087 for(int i = 0; i < _ReadInputnumw; i++)
0088 {
0089 _Farray[i] = new double[_ReadInputnumy];
0090 }
0091 }
0092 if(!_Farray1)
0093 {
0094 _Farray1 = new double*[_ReadInputnumw];
0095 for(int i = 0; i < _ReadInputnumw; i++)
0096 {
0097 _Farray1[i] = new double[_ReadInputnumy];
0098 }
0099 }
0100 if(!_Farray2)
0101 {
0102 _Farray2 = new double*[_ReadInputnumw];
0103 for(int i = 0; i < _ReadInputnumw; i++)
0104 {
0105 _Farray2[i] = new double[_ReadInputnumy];
0106 }
0107 }
0108
0109 double dummy[17];
0110
0111
0112 std::string wyFileName;
0113 wyFileName = _baseFileName +".txt";
0114
0115
0116
0117 double fpart =0.;
0118 double fptsum=0.;
0119 ifstream wylumfile;
0120
0121 _f_max=0.0;
0122 _f_max1=0.0;
0123 _f_max2=0.0;
0124
0125 wylumfile.open(wyFileName.c_str());
0126
0127 for(int i=0;i < 17;i++){
0128 wylumfile >> dummy[i];
0129 }
0130 int A_1 = dummy[1];
0131 int A_2 = dummy[3];
0132
0133 for(int i=0;i<_ReadInputnumw;i++){
0134 wylumfile >> _Warray[i];
0135 }
0136 for(int i=0;i<_ReadInputnumy;i++){
0137 wylumfile >> _Yarray[i];
0138 }
0139
0140 if( (_ReadInputgg_or_gP == 1) || (A_2 == 1 && A_1 != 1) || (A_1 ==1 && A_2 != 1) ){
0141 for(int i=0;i<_ReadInputnumw;i++){
0142 for(int j=0;j<_ReadInputnumy;j++){
0143 wylumfile >> _Farray[i][j];
0144 if( _Farray[i][j] > _f_max ) _f_max=_Farray[i][j];
0145 }
0146 }
0147
0148 for(int i=0;i<_ReadInputnumw;i++){
0149 for(int j=0;j<_ReadInputnumy;j++){
0150 _Farray[i][j]=_Farray[i][j]/_f_max;
0151 }
0152 }
0153 } else {
0154 for(int i=0;i<_ReadInputnumw;i++){
0155 for(int j=0;j<_ReadInputnumy;j++){
0156 wylumfile >> _Farray1[i][j];
0157 }
0158 }
0159 for(int i=0;i<_ReadInputnumw;i++){
0160 for(int j=0;j<_ReadInputnumy;j++){
0161 wylumfile >> _Farray2[i][j];
0162 if( _Farray1[i][j] + _Farray2[i][j] > _f_max ) _f_max=(_Farray1[i][j] + _Farray2[i][j]);
0163 }
0164 }
0165
0166 for(int i=0;i<_ReadInputnumw;i++){
0167 for(int j=0;j<_ReadInputnumy;j++){
0168 _Farray1[i][j]=_Farray1[i][j]/_f_max;
0169 _Farray2[i][j]=_Farray2[i][j]/_f_max;
0170 _Farray[i][j]=_Farray1[i][j]+_Farray2[i][j];
0171 }
0172 }
0173 }
0174 wylumfile >> _bwnormsave;
0175
0176 if (_ReadInputgg_or_gP != 1 && _ReadInputinterferencemode != 0) {
0177
0178 double **finterm = new double*[starlightLimits::MAXWBINS];
0179 for (int i = 0; i < starlightLimits::MAXWBINS; i++) finterm[i] = new double[starlightLimits::MAXYBINS];
0180 for (int i=0;i<_ReadInputnumy;i++) {
0181
0182
0183 fptsum=0.;
0184 for (int j=0;j<_ReadInputNPT;j++) {
0185 wylumfile >> fpart;
0186 finterm[i][j] = fpart;
0187 _fptarray[i][j]=0.;
0188 fptsum=fptsum+fpart;
0189 }
0190
0191 _fptarray[i][0]=finterm[i][0]/fptsum;
0192 for (int j=1;j<_ReadInputNPT;j++) {
0193 for (int k=0;k<j;k++) {
0194 _fptarray[i][j]=_fptarray[i][j]+finterm[i][k];
0195 }
0196 _fptarray[i][j]=_fptarray[i][j]/fptsum;
0197 }
0198 }
0199 delete [] finterm;
0200
0201 }
0202 wylumfile.close();
0203 return;
0204 }
0205
0206
0207
0208 void readLuminosity::e_read()
0209 {
0210 if(!_Warray) _Warray = new double[_ReadInputnumw];
0211 if(!_Yarray) _Yarray = new double[_ReadInputnumy];
0212 if(!_BWarray) _BWarray = new double[_ReadInputnumw];
0213 if(!_f_WYmax)
0214 {
0215 _f_WYarray = new double*[_ReadInputnumega];
0216 for(int i = 0; i < _ReadInputnumega; i++)
0217 {
0218 _f_WYarray[i] = new double[_ReadInputnumQ2];
0219 }
0220 }
0221 if(!_g_Earray)
0222 {
0223 _g_Earray = new double*[_ReadInputnumega];
0224 for(int i = 0; i < _ReadInputnumega; i++)
0225 {
0226 _g_Earray[i] = new double[_ReadInputnumQ2];
0227 }
0228 }
0229
0230 double dummy[13];
0231
0232
0233 std::string wyFileName;
0234 wyFileName = _baseFileName +".txt";
0235
0236
0237
0238 ifstream wylumfile;
0239
0240 _f_WYmax=0.0;
0241 _g_Emax=0.0;
0242
0243 wylumfile.open(wyFileName.c_str());
0244
0245 for(int i=0;i < 13;i++){
0246 wylumfile >> dummy[i];
0247 }
0248 int A_1 = dummy[1];
0249 int A_2 = dummy[3];
0250
0251 for(int i=0;i<_ReadInputnumw;i++){
0252 wylumfile >> _Warray[i];
0253 }
0254
0255 for(int i=0;i<_ReadInputnumy;i++){
0256 wylumfile >> _Yarray[i];
0257 }
0258
0259
0260 double bw_max = 0 ;
0261 for(int i=0;i<_ReadInputnumw;i++){
0262 wylumfile >> _BWarray[i];
0263 if( _BWarray[i] > bw_max )
0264 bw_max = _BWarray[i];
0265 }
0266 for(int i=0;i<_ReadInputnumw;i++){
0267 _BWarray[i] = _BWarray[i]/bw_max;
0268 }
0269
0270
0271 if( (A_2 == 0 && A_1 >= 1) || (A_1 ==0 && A_2 >= 1) ){
0272 for(int i=0;i<_ReadInputnumega;i++){
0273 for(int j=0;j<_ReadInputnumQ2;j++){
0274 wylumfile >> _f_WYarray[i][j];
0275 if( _f_WYarray[i][j] > _f_WYmax ) _f_WYmax=_f_WYarray[i][j];
0276
0277
0278
0279 }
0280 }
0281
0282 for(int i=0;i<_ReadInputnumega;i++){
0283 for(int j=0;j<_ReadInputnumQ2;j++){
0284 _f_WYarray[i][j] = _f_WYarray[i][j]/( _f_WYmax );
0285
0286 }
0287 }
0288 }
0289 wylumfile >> _bwnormsave;
0290
0291 wylumfile.close();
0292 cout<<"Done reading wylumi file"<<endl;
0293 std::string EQ2FileName;
0294 EQ2FileName = "e_"+_baseFileName+".txt";
0295 ifstream EQlumfile;
0296
0297 EQlumfile.open(EQ2FileName.c_str());
0298 int n_steps;
0299 EQlumfile >> n_steps;
0300 double integrated_max = 0 ;
0301
0302 _g_EQ2array = new vector< std::pair< double, vector<double> > > ();
0303 while( !EQlumfile.eof() ){
0304 _g_EQ2max = 0 ;
0305 double integral;
0306 std::vector<double> p;
0307 for( int iQ2 = 0 ; iQ2 < _ReadInputnumQ2+3 ; ++iQ2){
0308 if(iQ2 == 0 ){
0309 EQlumfile >> integral;
0310 if( integral > integrated_max)
0311 integrated_max = integral;
0312 }
0313 else{
0314 double temp;
0315 EQlumfile >> temp;
0316 p.push_back(temp);
0317
0318 }
0319 if( iQ2 > 2 && p[iQ2-1] > _g_EQ2max){
0320 _g_EQ2max = p[iQ2-1];
0321 }
0322 }
0323 for( unsigned int iQ2=2; iQ2 < p.size(); ++iQ2)
0324 p[iQ2] = p[iQ2]/_g_EQ2max;
0325
0326 _g_EQ2array->push_back(std::pair< double, std::vector<double> >(integral,p));
0327 }
0328 for( std::vector< std::pair<double, std::vector<double> > >::iterator it =_g_EQ2array->begin() ; it != _g_EQ2array->end(); ++it){
0329
0330 it->first = it->first/integrated_max;
0331
0332 }
0333
0334
0335
0336
0337
0338
0339
0340
0341
0342 EQlumfile.close();
0343 return;
0344
0345 }
0346