estarlight/src/readinluminosity.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 265                         $: revision of last commit
0024 // $Author:: butter                   $: author of last commit
0025 // $Date:: 2016-06-06 21:37:26 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //    Added 18->19 for reading in the luminosity table
0029 //    Incoherent factor added to table --Joey
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)/*, _g_EQ2array(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   // For eSTARLIGHT
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]; //number of lines used to read in input parameters saved to lookup table[slight.txt].
0110
0111
0112   std::string wyFileName;
0113   wyFileName = _baseFileName +".txt";
0114
0115 //  cout << "wyFileName being read in" << wyFileName << endl;
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     //Normalize farray
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     //Normalize farray, farray1, farray2
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         // only numy/2 y bins here, from 0 (not -ymax) to ymax
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             //fmax=0;
0182             //we want to convert _fptarray to an integral array where fpt(i,j) is near 0, and fpt(j,NPT) ~1. This will facilitate a simple table loookup
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             //convert array to integral
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]; //number of lines used to read in input parameters saved to lookup table[slight.txt].
0231
0232
0233   std::string wyFileName;
0234   wyFileName = _baseFileName +".txt";
0235
0236 //  cout << "wyFileName being read in" << wyFileName << endl;
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   //No longer needed, can be removed later
0255   for(int i=0;i<_ReadInputnumy;i++){
0256     wylumfile >> _Yarray[i];
0257   }
0258
0259   //New table saving the Breit-Wigner function for form factor
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     //wylumfile >> _g_Earray[i][j];
0278     //if( _g_Earray[i][j] > _g_Emax ) _g_Emax = _g_Earray[i][j];
0279       }
0280     }
0281     //Normalize f_WY array, g does not need to be normalized, it is used for normalization
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     //_g_Earray[i][j] = _g_Earray[i][j]/_g_Emax;
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   //  _g_EQ2array = new map<string,std::vector<double> >();
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     //cout<<p[iQ2-1]<<endl;
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     //cout<<it->first<<endl;
0330     it->first = it->first/integrated_max;
0331     //cout<<"Now"<<it->first<<endl;
0332   }
0333
0334   //normalize
0335   /*for(  std::map<string,std::vector<double> >::iterator it =_g_EQ2array->begin() ; it != _g_EQ2array->end(); ++it){
0336     std::string key = it->first;
0337     if(key=="")
0338       continue;
0339     for(unsigned int iQ2=2; iQ2 < it->second.size(); ++iQ2)
0340       it->second[iQ2] = it->second[iQ2]/_g_EQ2max;
0341       }*/
0342   EQlumfile.close();
0343   return;
0344
0345 }
0346