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