Completing the setup for TTherminator model.
[u/mrichter/AliRoot.git] / TTherminator / Therminator / Hypersurface.cxx
1 #include "Hypersurface.h"
2
3 Hypersurface::Hypersurface(const char *dirname) {
4 // Read FOHSI.txt file and init parameters
5   char *dirfull;
6   dirfull = (char *) malloc(sizeof(char)*strlen(dirname) + 5);
7   sprintf(dirfull, "%s", dirname);
8   if (strlen(dirfull) > 0) 
9     if (dirfull[strlen(dirfull)-1] != '/')
10       strcat(dirfull, "/");
11   FName = (char *) malloc(sizeof(char)* strlen(dirfull)+ 50);
12   sprintf(FName,"%sFOHSI.txt", dirfull);
13   HSFile = new ifstream;
14   HSFile->open(FName);
15   if(HSFile->is_open()) {
16     HSFile->seekg (0, ios::beg);
17 // phi
18     HSFile->getline(buff,100);
19     Np=atoi(buff);
20     HSFile->getline(buff,100);
21     ip=atof(buff);
22     HSFile->getline(buff,100);
23     fp=atof(buff);
24 // zeta
25     HSFile->getline(buff,100);
26     Nz=atoi(buff);
27     HSFile->getline(buff,100);
28     iz=atof(buff);
29     HSFile->getline(buff,100);
30     fz=atof(buff);
31 // freeze-out temperature
32     HSFile->getline(buff,100);
33     TFO=atof(buff);
34 // hydro initial time
35     HSFile->getline(buff,100);
36     tau0=atof(buff);
37
38     HSFile->close();
39     dp = (fp-ip)/(Np-1);
40     dz = (fz-iz)/(Nz-1);
41
42 // create 2D array
43     aArr   = new double* [Np];
44     vArr   = new double* [Np];
45     dArr   = new double* [Np];
46     DpdArr = new double* [Np];
47     DzdArr = new double* [Np];
48     for(i=0;i<Np;i++) {
49       aArr[i]   = new double [Nz];
50       vArr[i]   = new double [Nz];
51       dArr[i]   = new double [Nz];
52       DpdArr[i] = new double [Nz];
53       DzdArr[i] = new double [Nz];
54     }
55   } else {
56     Np = 0; ip = 0.0; fp = 0.0; dp = 1.0;
57     Nz = 0; iz = 0.0; fz = 0.0; dz = 1.0;
58     TFO = 0.0;
59     aArr   = NULL;
60     vArr   = NULL;
61     dArr   = NULL;
62     DpdArr = NULL;
63     DzdArr = NULL;
64   }
65   delete HSFile;
66
67   // Read FOHSa.txt file and fill array
68   sprintf(FName,"%sFOHSa.txt", dirfull);
69   //  FName = "FOHSa.txt";
70   HSFile = new ifstream;
71   HSFile->open(FName);
72   if(HSFile->is_open()) {
73     HSFile->seekg (0, ios::beg);
74     for(i=0;i<Np;i++) {
75       for(j=0;j<Nz;j++) {
76         HSFile->getline(buff,100);
77         aArr[i][j]=atof(buff);
78       }
79     }
80     HSFile->close();
81   }
82   delete HSFile;
83
84   // Read FOHSv.txt file and fill array
85   sprintf(FName,"%sFOHSv.txt", dirfull);
86   //  FName = "FOHSv.txt";
87   HSFile = new ifstream;
88   HSFile->open(FName);
89   if(HSFile->is_open()) {
90     HSFile->seekg (0, ios::beg);
91     for(i=0;i<Np;i++) {
92       for(j=0;j<Nz;j++) {
93         HSFile->getline(buff,100);
94         vArr[i][j]=atof(buff);
95       }
96     }
97     HSFile->close();
98   }
99   delete HSFile;
100
101   // Read FOHSd.txt file and fill array
102   sprintf(FName,"%sFOHSd.txt", dirfull);
103   //  FName = "FOHSd.txt";
104   HSFile = new ifstream;
105   HSFile->open(FName);
106   if(HSFile->is_open()) {
107     HSFile->seekg (0, ios::beg);
108     for(i=0;i<Np;i++) {
109       for(j=0;j<Nz;j++) {
110         HSFile->getline(buff,100);
111         dArr[i][j]=atof(buff);
112       }
113     }
114     HSFile->close();
115   }
116   delete HSFile;
117
118   // Read FOHSDpd.txt file and fill array
119   sprintf(FName,"%sFOHSDpd.txt", dirfull);
120   //  FName = "FOHSDpd.txt";
121   HSFile = new ifstream;
122   HSFile->open(FName);
123   if(HSFile->is_open()) {
124     HSFile->seekg (0, ios::beg);
125     for(i=0;i<Np;i++) {
126       for(j=0;j<Nz;j++) {
127         HSFile->getline(buff,100);
128         DpdArr[i][j]=atof(buff);
129       }
130     }
131     HSFile->close();
132   }
133   delete HSFile;
134
135   // Read FOHSDzd.txt file and fill array
136   sprintf(FName,"%sFOHSDzd.txt", dirfull);
137   //  FName = "FOHSDzd.txt";
138   HSFile = new ifstream;
139   HSFile->open(FName);
140   if(HSFile->is_open()) {
141     HSFile->seekg (0, ios::beg);
142     for(i=0;i<Np;i++) {
143       for(j=0;j<Nz;j++) {
144         HSFile->getline(buff,100);
145         DzdArr[i][j]=atof(buff);
146       }
147     }
148     HSFile->close();
149   }
150   delete HSFile;
151   free (FName);
152   free (dirfull);
153 }
154
155 Hypersurface::Hypersurface(void) {
156   Hypersurface::Hypersurface("./");
157 }
158
159
160 Hypersurface::~Hypersurface(void) {
161   for(i=0;i<Np;i++) {
162     delete[] aArr[i];
163     delete[] vArr[i];
164     delete[] dArr[i];
165     delete[] DpdArr[i];
166     delete[] DzdArr[i];
167   }
168   delete[] aArr;
169   delete[] vArr;
170   delete[] dArr;
171   delete[] DpdArr;
172   delete[] DzdArr;
173 }
174
175 double Hypersurface::fahs(double p, double z) {
176   i=int(p/dp);
177   j=int(z/dz);
178   if(i>=Np-1) i=Np-2;
179   if(j>=Nz-1) j=Nz-2;
180   return (aArr[i][j]   * (i+1-p/dp) + aArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
181          (aArr[i][j+1] * (i+1-p/dp) + aArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
182 }
183
184 double Hypersurface::fvhs(double p, double z) {
185   i=int(p/dp);
186   j=int(z/dz);
187   if(i>=Np-1) i=Np-2;
188   if(j>=Nz-1) j=Nz-2;
189   return (vArr[i][j]   * (i+1-p/dp) + vArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
190          (vArr[i][j+1] * (i+1-p/dp) + vArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
191 }
192
193 double Hypersurface::fdhs(double p, double z) {
194   i=int(p/dp);
195   j=int(z/dz);
196   if(i>=Np-1) i=Np-2;
197   if(j>=Nz-1) j=Nz-2;
198   return (dArr[i][j]   * (i+1-p/dp) + dArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
199          (dArr[i][j+1] * (i+1-p/dp) + dArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
200 }
201
202 double Hypersurface::fDpdhs(double p, double z) {
203   i=int(p/dp);
204   j=int(z/dz);
205   if(i>=Np-1) i=Np-2;
206   if(j>=Nz-1) j=Nz-2;
207   return (DpdArr[i][j]   * (i+1-p/dp) + DpdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
208          (DpdArr[i][j+1] * (i+1-p/dp) + DpdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
209 }
210
211 double Hypersurface::fDzdhs(double p, double z) {
212   i=int(p/dp);
213   j=int(z/dz);
214   if(i>=Np-1) i=Np-2;
215   if(j>=Nz-1) j=Nz-2;
216   return (DzdArr[i][j]   * (i+1-p/dp) + DzdArr[i+1][j]   * (p/dp-i)) * (j+1-z/dz) +
217          (DzdArr[i][j+1] * (i+1-p/dp) + DzdArr[i+1][j+1] * (p/dp-i)) * (z/dz-j);
218 }