STARLIGHT/starlight/src/readinluminosity.cpp

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 //    Copyright 2010
   4 //
   5 //    This file is part of starlight.
   6 //
   7 //    starlight is free software: you can redistribute it and/or modify
   8 //    it under the terms of the GNU General Public License as published by
   9 //    the Free Software Foundation, either version 3 of the License, or
  10 //    (at your option) any later version.
  11 //
  12 //    starlight is distributed in the hope that it will be useful,
  13 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 //    GNU General Public License for more details.
  16 //
  17 //    You should have received a copy of the GNU General Public License
  18 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
  19 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev:: 193                         $: revision of last commit
  24 // $Author:: jnystrand                $: author of last commit
  25 // $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
  26 //
  27 // Description:
  28 //    Added 18->19 for reading in the luminosity table
  29 //    Incoherent factor added to table --Joey
  30 //
  31 //
  32 //
  33 ///////////////////////////////////////////////////////////////////////////
  34
  35
  36 #include <iostream>
  37 #include <fstream>
  38
  39 #include "readinluminosity.h"
  40 #include "starlightconstants.h"
  41 #include "inputParameters.h"
  42
  43 using namespace std;
  44
  45
  46 //______________________________________________________________________________
  47 readLuminosity::readLuminosity()//:inputread(input)
  48   : _Warray(0), _Yarray(0), _Farray(0), _Farray1(0), _Farray2(0)
  49 {
  50   //storing inputparameters into protected variables for the object to use them
  51   _ReadInputNPT=inputParametersInstance.nmbPtBinsInterference();
  52   _ReadInputnumy=inputParametersInstance.nmbRapidityBins();
  53   _ReadInputnumw=inputParametersInstance.nmbWBins();
  54   _ReadInputgg_or_gP=inputParametersInstance.productionMode();
  55   _ReadInputinterferencemode=inputParametersInstance.interferenceEnabled();
  56
  57 }
  58
  59
  60 //______________________________________________________________________________
  61 readLuminosity::~readLuminosity()
  62 {
  63   if(_Warray) delete [] _Warray;
  64   if(_Yarray) delete [] _Yarray;
  65   if(_Farray) delete [] _Farray;
  66   if(_Farray1) delete [] _Farray1;
  67   if(_Farray2) delete [] _Farray2;
  68 }
  69
  70
  71 //______________________________________________________________________________
  72 void readLuminosity::read()
  73 {
  74
  75   if(!_Warray) _Warray = new double[_ReadInputnumw];
  76   if(!_Yarray) _Yarray = new double[_ReadInputnumy];
  77   if(!_Farray)
  78   {
  79     _Farray = new double*[_ReadInputnumw];
  80     for(int i = 0; i < _ReadInputnumw; i++)
  81     {
  82       _Farray[i] = new double[_ReadInputnumy];
  83     }
  84   }
  85   if(!_Farray1)
  86   {
  87     _Farray1 = new double*[_ReadInputnumw];
  88     for(int i = 0; i < _ReadInputnumw; i++)
  89     {
  90       _Farray1[i] = new double[_ReadInputnumy];
  91     }
  92   }
  93   if(!_Farray2)
  94   {
  95     _Farray2 = new double*[_ReadInputnumw];
  96     for(int i = 0; i < _ReadInputnumw; i++)
  97     {
  98       _Farray2[i] = new double[_ReadInputnumy];
  99     }
 100   }
 101
 102   double dummy[19]; //14//18
 103
 104   //  double (*finterm)[starlightLimits::MAXWBINS]=new double[starlightLimits::MAXWBINS][starlightLimits::MAXYBINS];
 105   //decreased from 1000*1000; too big! causes fault!
 106
 107   double fpart =0.;
 108   double fptsum=0.;
 109   ifstream wylumfile;
 110
 111   _f_max=0.0;
 112   _f_max1=0.0;
 113   _f_max2=0.0;
 114
 115   wylumfile.open("slight.txt");
 116   for(int i=0;i < 19;i++){ // was 14; this is to account for sergei's additional parameters ie d-Au//was19
 117     wylumfile >> dummy[i];
 118   }
 119   int A_1 = int(dummy[1]);
 120   int A_2 = int(dummy[3]);
 121
 122   for(int i=0;i<_ReadInputnumw;i++){
 123     wylumfile >> _Warray[i];
 124   }
 125   for(int i=0;i<_ReadInputnumy;i++){
 126     wylumfile >> _Yarray[i];
 127   }
 128
 129   if( (_ReadInputgg_or_gP == 1) || (A_2 == 1 && A_1 != 1) || (A_1 ==1 && A_2 != 1) ){
 130     for(int i=0;i<_ReadInputnumw;i++){
 131       for(int j=0;j<_ReadInputnumy;j++){
 132         wylumfile >> _Farray[i][j];
 133         if( _Farray[i][j] > _f_max ) _f_max=_Farray[i][j];
 134       }
 135     }
 136     //Normalize farray
 137     for(int i=0;i<_ReadInputnumw;i++){
 138       for(int j=0;j<_ReadInputnumy;j++){
 139         _Farray[i][j]=_Farray[i][j]/_f_max;
 140       }
 141     }
 142   } else {
 143     for(int i=0;i<_ReadInputnumw;i++){
 144       for(int j=0;j<_ReadInputnumy;j++){
 145         wylumfile >> _Farray1[i][j];
 146         //        if( _Farray1[i][j] > _f_max ) _f_max=_Farray1[i][j];
 147       }
 148     }
 149     for(int i=0;i<_ReadInputnumw;i++){
 150       for(int j=0;j<_ReadInputnumy;j++){
 151         wylumfile >> _Farray2[i][j];
 152         if( _Farray1[i][j] + _Farray2[i][j] > _f_max ) _f_max=(_Farray1[i][j] + _Farray2[i][j]);
 153       }
 154     }
 155     //Normalize farray, farray1, farray2
 156     for(int i=0;i<_ReadInputnumw;i++){
 157       for(int j=0;j<_ReadInputnumy;j++){
 158         _Farray1[i][j]=_Farray1[i][j]/_f_max;
 159         _Farray2[i][j]=_Farray2[i][j]/_f_max;
 160         _Farray[i][j]=_Farray1[i][j]+_Farray2[i][j];
 161       }
 162     }
 163   }
 164
 165   if (_ReadInputgg_or_gP != 1 && _ReadInputinterferencemode != 0) {
 166         // only numy/2 y bins here, from 0 (not -ymax) to ymax
 167         double **finterm  = new double*[starlightLimits::MAXWBINS];
 168         for (int i = 0; i < starlightLimits::MAXWBINS; i++) finterm[i] = new double[starlightLimits::MAXYBINS];
 169         for (int i=0;i<_ReadInputnumy/2;i++) {
 170             //fmax=0;
 171             //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
 172             fptsum=0.;
 173             for (int j=0;j<_ReadInputNPT;j++) {
 174                 wylumfile >> fpart;
 175                 finterm[i][j] = fpart;
 176                 _fptarray[i][j]=0.;
 177                 fptsum=fptsum+fpart;
 178             }
 179             //convert array to integral
 180             _fptarray[i][0]=finterm[i][0]/fptsum;
 181             for (int j=1;j<_ReadInputNPT;j++) {
 182                 for (int k=0;k<=j;k++) {
 183                     _fptarray[i][j]=_fptarray[i][j]+finterm[i][k];
 184                 }
 185                 _fptarray[i][j]=_fptarray[i][j]/fptsum;
 186             }
 187         }
 188         delete [] finterm;
 189
 190     }
 191   wylumfile >> _bwnormsave;
 192   wylumfile.close();
 193   //8delete[] finterm;
 194   return;
 195 }