]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - THbtp/THBTprocessor.cxx
Transition to NewIO
[u/mrichter/AliRoot.git] / THbtp / THBTprocessor.cxx
... / ...
CommitLineData
1#include <TGenerator.h>
2#include "THBTprocessor.h"
3#include <TParticle.h>
4#include <TMath.h>
5
6#include <Riostream.h>
7#ifndef WIN32
8# define hbtprocessor hbtprocessor_
9# define ctest ctest_
10# define type_of_call
11#else
12# define hbtprocessor HBTPROCESSOR
13# define ctest CTEST
14# define type_of_call _stdcall
15#endif
16
17
18ClassImp(THBTprocessor)
19
20extern "C" void type_of_call hbtprocessor();
21extern "C" void type_of_call ctest();
22
23//Wrapper class to HBT processor fortran code.
24//For more information see AliGenHBTprocessor class
25//HBT processor is written by Lanny Ray
26//Wrapper class written by Piotr Krzysztof Skowronski (Piotr.Skowronski@cern.ch)
27//
28
29/*************************************************/
30THBTprocessor::THBTprocessor()// it is better not to intialize it:TGenerator("THBTprocessor","THBTprocessor")
31 //because it allocates memmory for TObjArray::fParticles which is not used in our case
32 //and we are using TClonesArray in import paerticles
33 {
34 //constructor
35 PARAMETERS.ALICE = 1; //flag that we are working in AliRoot (0==STAR stand alone mode)
36 PARAMETERS.pi = TMath::Pi();//3.141592654;
37 PARAMETERS.rad = 180.0/TMath::Pi();
38 PARAMETERS.hbc = 0.19732891;
39
40 fParticles = 0;
41 Initialize(); //Enforce initialization (cleaning all commons)
42 }
43
44void THBTprocessor::GenerateEvent()
45{
46//Starts processing
47
48 cout<<endl<<"Entering Fortran"<<endl;
49 ctest();
50 hbtprocessor();
51 cout<<endl<<"Exited Fortran"<<endl;
52 if(PARAMETERS.errorcode)
53 {
54 Fatal("THBTprocessor::GenerateEvent()","HBT Processor (fortran part) finished with errors\n \
55 See hbt_simulation.out file for more detailed information");
56 }
57 else
58 cout<<endl<<"GOOD ! HBT Processor finished without errors"<<endl;
59 cout<<"PARAMETERS.errorcode= "<<PARAMETERS.errorcode<<endl;
60
61}
62/*****************************************************************************************/
63void THBTprocessor::Initialize()
64{
65 //IT RESETS ALL THE PREVIOUS SETTINGS
66 //Call this method to set default values in PARAMETERS & MESH
67 //and zero other common block
68
69 if(gDebug) cout<<"\nHBT PROCESSOR::Initialize() Setting Default valuses in all COMMON BLOCKS"<<endl;
70
71 PARAMETERS.ref_control = 2;
72 PARAMETERS.switch_1d = 3;
73 PARAMETERS.switch_3d = 1;
74 PARAMETERS.switch_type = 3;
75 PARAMETERS.switch_coherence = 0;
76 PARAMETERS.switch_coulomb = 3;
77 PARAMETERS.switch_fermi_bose= 1;
78 PARAMETERS.trk_accep = 1.0;
79 PARAMETERS.print_full = 1;
80 PARAMETERS.print_sector_data= 1;
81
82 PARAMETERS.n_pid_types = 2;
83 PARAMETERS.pid[0] = 8;
84 PARAMETERS.pid[1] = 9;
85 PARAMETERS.deltap = 0.1;
86 PARAMETERS.maxit = 50;
87 PARAMETERS.delchi = 1.0;
88 PARAMETERS.irand = 76564;
89 PARAMETERS.lambda = 0.6;
90 PARAMETERS.R_1d = 7.0;
91 PARAMETERS.Rside = 6.0;
92
93 PARAMETERS.Rout = 7.0;
94 PARAMETERS.Rlong = 4.0;
95 PARAMETERS.Rperp = 6.0;
96 PARAMETERS.Rparallel = 4.0;
97 PARAMETERS.R0 = 4.0;
98 PARAMETERS.Q0 = 9.0;
99
100 PARAMETERS.n_part_1_trk = 0;
101 PARAMETERS.n_part_2_trk = 0;
102 PARAMETERS.n_part_tot_trk = 0;
103 PARAMETERS.n_part_used_1_trk = 0;
104
105 PARAMETERS.n_part_used_2_trk = 0;
106 PARAMETERS.n_part_1_trk2 = 0;
107 PARAMETERS.n_part_2_trk2 = 0;
108 PARAMETERS.n_part_tot_trk2 = 0;
109 PARAMETERS.n_part_used_1_trk2 = 0;
110 PARAMETERS.n_part_used_2_trk2 = 0;
111 PARAMETERS.n_part_used_1_ref = 0;
112 PARAMETERS.n_part_used_2_ref = 0;
113 PARAMETERS.n_part_used_1_inc = 0;
114 PARAMETERS.n_part_used_2_inc = 0;
115
116 PARAMETERS.num_pairs_like = 0;
117 PARAMETERS.num_pairs_unlike = 0;
118 PARAMETERS.num_pairs_like_ref = 0;
119 PARAMETERS.num_pairs_like_inc = 0;
120 PARAMETERS.num_pairs_unlike_inc = 0;
121 PARAMETERS.event_line_counter = 0;
122 PARAMETERS.file10_line_counter = 0;
123
124 PARAMETERS.chisq_wt_like_1d = 1.0;
125 PARAMETERS.chisq_wt_unlike_1d = 1.0;
126 PARAMETERS.chisq_wt_like_3d_fine = 1.0;
127
128 PARAMETERS.chisq_wt_unlike_3d_fine = 1.0;
129 PARAMETERS.chisq_wt_like_3d_coarse = 1.0;
130 PARAMETERS.chisq_wt_unlike_3d_coarse= 1.0;
131 PARAMETERS.chisq_wt_hist1_1 = 1.0;
132 PARAMETERS.chisq_wt_hist1_2 = 1.0; // /////internal comment 25 fields
133
134/*********************************************/
135
136
137 MESH.n_pt_bins = 50; //OK
138 MESH.pt_min = 0.1; //Pt in GeV/c //OK
139 MESH.pt_max = 0.98; //Pt in GeV/c
140 MESH.n_phi_bins = 50; //OK
141 MESH.phi_min = 0.0; //OK
142 MESH.phi_max = 360.0; //OK
143 MESH.n_eta_bins = 50; //OK
144 MESH.eta_min =-1.5; //OK
145 MESH.eta_max = 1.5; //OK
146 MESH.n_1d_fine = 10; //OK
147 MESH.binsize_1d_fine = 0.01; //ok
148 MESH.n_1d_coarse = 2; //O
149 MESH.binsize_1d_coarse= 0.05; //ok
150 MESH.n_3d_fine = 8; //OK
151 MESH.binsize_3d_fine = 0.01; //ok
152 MESH.n_3d_coarse = 2; //OK
153 MESH.binsize_3d_coarse= 0.08; //ok
154 MESH.n_3d_fine_project= 3; //OK
155 MESH.n_px_bins = 20; //OK
156 MESH.px_min =-1.0; //OK
157 MESH.px_max = 1.0; //OK
158 MESH.n_py_bins = 20; //OK
159 MESH.py_min =-1.0; //OK
160 MESH.py_max = 1.0; //OK
161 MESH.n_pz_bins = 70; //OK
162 MESH.pz_min =-3.6; //OK
163 MESH.pz_max = 3.6; //OK
164
165/*********************************************/
166
167 Int_t i; //loop variable
168
169 for (i =0; i<TRK_MAXLEN;i++)
170 {
171 TRACK1.trk_id[i] = 0;
172 TRACK1.trk_px_sec[i] = 0;
173 TRACK1.trk_py_sec[i] = 0;
174 TRACK1.trk_pz_sec[i] = 0;
175 TRACK1.trk_sector[i] = 0;
176 TRACK1.trk_flag[i] = 0;
177 TRACK1.trk_out_flag[i] = 0;
178 TRACK1.trk_merge_flag[i] = 0;
179 TRACK1.trk_ge_pid[i] = 0;
180 TRACK1.trk_start_vertex[i] = 0;
181 TRACK1.trk_stop_vertex[i] = 0;
182 TRACK1.trk_event_line[i] = 0;
183 TRACK1.trk_px[i] = 0;
184 TRACK1.trk_py[i] = 0;
185 TRACK1.trk_pz[i] = 0;
186 TRACK1.trk_E[i] = 0;
187 TRACK1.trk_pt[i] = 0;
188 TRACK1.trk_phi[i] = 0;
189 TRACK1.trk_eta[i] = 0;
190 }
191
192/*********************************************/
193
194 for (i =0; i<TRK2_MAXLEN;i++)
195 {
196 TRACK2.trk2_id[i] = 0;
197 TRACK2.trk2_px_sec[i] = 0;
198 TRACK2.trk2_py_sec[i] = 0;
199 TRACK2.trk2_pz_sec[i] = 0;
200 TRACK2.trk2_sector[i] = 0;
201 TRACK2.trk2_flag[i] = 0;
202 TRACK2.trk2_out_flag[i] = 0;
203 TRACK2.trk2_merge_flag[i] = 0;
204 TRACK2.trk2_ge_pid[i] = 0;
205 TRACK2.trk2_start_vertex[i] = 0;
206 TRACK2.trk2_stop_vertex[i] = 0;
207 TRACK2.trk2_event_line[i] = 0;
208 TRACK2.trk2_px[i] = 0;
209 TRACK2.trk2_py[i] = 0;
210 TRACK2.trk2_pz[i] = 0;
211 TRACK2.trk2_E[i] = 0;
212 TRACK2.trk2_pt[i] = 0;
213 TRACK2.trk2_phi[i] = 0;
214 TRACK2.trk2_eta[i] = 0;
215 }
216
217/*********************************************/
218
219 for (i =0; i<SEC_MAXLEN;i++)
220 {
221 SEC_TRK_MAP.stm_sec_id [i] = 0;
222 SEC_TRK_MAP.stm_n_trk_sec[i] = 0;
223 SEC_TRK_MAP.stm_flag[i] = 0;
224
225 for (Int_t j=0; j<MAX_TRK_SEC;j++)
226 SEC_TRK_MAP.stm_track_id[i][j] = 0;
227 }
228
229/*********************************************/
230
231 for (i =0; i<SEC_MAXLEN2;i++)
232 {
233 SEC_TRK_MAP2.stm_sec_id2[i] = 0;
234 SEC_TRK_MAP2.stm_n_trk_sec2[i] = 0;
235 SEC_TRK_MAP2.stm_flag2[i] = 0;
236
237 for (Int_t j=0; j<MAX_TRK_SEC;j++)
238 SEC_TRK_MAP2.stm_track_id2[i][j] = 0;
239 }
240
241/*********************************************/
242
243 for (i =0; i<PART_MAXLEN;i++)
244 {
245 PARTICLE.part_id[i] = 0;
246 PARTICLE.part_charge[i] = 0;
247 PARTICLE.part_mass[i] = 0;
248 PARTICLE.part_lifetime[i] = 0;
249 }
250
251
252/*********************************************/
253 for (i =0; i<MAX_C2_1D;i++)
254 {
255 CORRELATIONS.c2mod_like_1d[i] = 0;
256 CORRELATIONS.c2mod_unlike_1d[i] = 0;
257 CORRELATIONS.c2fit_like_1d[i] = 0;
258 CORRELATIONS.c2fit_unlike_1d[i] = 0;
259 CORRELATIONS.c2err_like_1d[i] = 0;
260 CORRELATIONS.c2err_unlike_1d[i] = 0;
261 }
262/*********************************************/
263 for (i =0; i<MAX_C2_3D;i++)
264 for (Int_t j =0; j<MAX_C2_3D;j++)
265 for (Int_t k =0; k<MAX_C2_3D;k++)
266 {
267 CORRELATIONS.c2mod_like_3d_fine[i][j][k] = 0.0;
268 CORRELATIONS.c2mod_unlike_3d_fine[i][j][k] = 0.0;
269 CORRELATIONS.c2mod_like_3d_coarse[i][j][k] = 0.0;
270 CORRELATIONS.c2mod_unlike_3d_coarse[i][j][k] = 0.0;
271 CORRELATIONS.c2fit_like_3d_fine[i][j][k] = 0.0;
272 CORRELATIONS.c2fit_unlike_3d_fine[i][j][k] = 0.0;
273 CORRELATIONS.c2fit_like_3d_coarse[i][j][k] = 0.0;
274 CORRELATIONS.c2fit_unlike_3d_coarse[i][j][k] = 0.0;
275 CORRELATIONS.c2err_like_3d_fine[i][j][k] = 0.0;
276 CORRELATIONS.c2err_unlike_3d_fine[i][j][k] = 0.0;
277 CORRELATIONS.c2err_like_3d_coarse[i][j][k] = 0.0;
278 CORRELATIONS.c2err_unlike_3d_coarse[i][j][k] = 0.0;
279 }
280/*********************************************/
281
282 EVENT_SUMMARY.niter_mean = 0.0;
283 EVENT_SUMMARY.niter_rms = 0.0;
284
285 EVENT_SUMMARY.npart1_mean = 0.0;
286 EVENT_SUMMARY.npart1_rms = 0.0;
287
288 EVENT_SUMMARY.npart2_mean = 0.0;
289 EVENT_SUMMARY.npart2_rms = 0.0;
290
291 EVENT_SUMMARY.npart_tot_mean = 0.0;
292 EVENT_SUMMARY.npart_tot_rms = 0.0;
293
294 EVENT_SUMMARY.nsec_flag_mean = 0.0;
295 EVENT_SUMMARY.nsec_flag_rms = 0.0;
296
297 EVENT_SUMMARY.frac_trks_out_mean = 0.0;
298 EVENT_SUMMARY.frac_trks_out_rms = 0.0;
299
300 EVENT_SUMMARY.frac_trks_flag_mean = 0.0;
301 EVENT_SUMMARY.frac_trks_flag_rms = 0.0;
302
303 EVENT_SUMMARY.chi_l1d_mean = 0.0;
304 EVENT_SUMMARY.chi_l1d_rms = 0.0;
305
306 EVENT_SUMMARY.chi_u1d_mean = 0.0;
307 EVENT_SUMMARY.chi_u1d_rms = 0.0;
308
309 for (i =0; i<MAX_EVENTS;i++)
310 {
311 EVENT_SUMMARY.n_part_used_1_store[i] = 0.0;
312 EVENT_SUMMARY.n_part_used_2_store[i] = 0.0;
313 EVENT_SUMMARY.n_part_tot_store[i] = 0.0;
314 EVENT_SUMMARY.num_sec_flagged_store[i] = 0.0;
315 EVENT_SUMMARY.frac_trks_out[i] = 0.0;
316 EVENT_SUMMARY.frac_trks_flag[i] = 0.0;
317 EVENT_SUMMARY.chisq_like_1d_store[i] = 0.0;
318 EVENT_SUMMARY.num_iter[i] = 0.0;
319 EVENT_SUMMARY.chisq_unlike_1d_store[i] = 0.0;
320 }
321/*********************************************/
322 for (i =0; i<MAX_C2_COUL;i++)
323 {
324 COULMB.c2_coul_like[i] = 0.0;
325 COULMB.c2_coul_unlike[i] = 0.0;
326 COULMB.q_coul[i] = 0.0;
327
328 }
329/*********************************************/
330 for (i =0; i<MAX_H_1D;i++)
331 {
332 HISTOGRAMS.hist1_pt_1[i] = 0;
333 HISTOGRAMS.hist1_phi_1[i] = 0;
334 HISTOGRAMS.hist1_eta_1[i] = 0;
335 HISTOGRAMS.hist1_pt_2[i] = 0;
336 HISTOGRAMS.hist1_phi_2[i] = 0;
337 HISTOGRAMS.hist1_eta_2[i] = 0;
338 HISTOGRAMS.htmp1_pt_1[i] = 0;
339 HISTOGRAMS.htmp1_phi_1[i] = 0;
340 HISTOGRAMS.htmp1_eta_1[i] = 0;
341 HISTOGRAMS.htmp1_pt_2[i] = 0;
342 HISTOGRAMS.htmp1_phi_2[i] = 0;
343 HISTOGRAMS.htmp1_eta_2[i] = 0;
344 HISTOGRAMS.href1_pt_1[i] = 0;
345 HISTOGRAMS.href1_phi_1[i] = 0;
346 HISTOGRAMS.href1_eta_1[i] = 0;
347 HISTOGRAMS.href1_pt_2[i] = 0;
348 HISTOGRAMS.href1_phi_2[i] = 0;
349 HISTOGRAMS.href1_eta_2[i] = 0;
350 HISTOGRAMS.hinc1_pt_1[i] = 0;
351 HISTOGRAMS.hinc1_phi_1[i] = 0;
352 HISTOGRAMS.hinc1_eta_1[i] = 0;
353 HISTOGRAMS.hinc1_pt_2[i] = 0;
354 HISTOGRAMS.hinc1_phi_2[i] = 0;
355 HISTOGRAMS.hinc1_eta_2[i] = 0;
356 HISTOGRAMS.hist_like_1d[i] = 0;
357 HISTOGRAMS.hist_unlike_1d[i] = 0;
358 HISTOGRAMS.htmp_like_1d[i] = 0;
359 HISTOGRAMS.htmp_unlike_1d[i] = 0;
360 HISTOGRAMS.href_like_1d[i] = 0;
361 HISTOGRAMS.href_unlike_1d[i] = 0;
362 HISTOGRAMS.hinc_like_1d[i] = 0;
363 HISTOGRAMS.hinc_unlike_1d[i] = 0;
364 }
365
366 for (i =0; i<MAX_H_3D;i++)
367 for (Int_t j =0; j<MAX_H_3D;j++)
368 for (Int_t k =0; k<MAX_H_3D;k++)
369 {
370 HISTOGRAMS.hist_like_3d_fine[i][j][k] = 0;
371 HISTOGRAMS.hist_unlike_3d_fine[i][j][k] = 0;
372 HISTOGRAMS.hist_like_3d_coarse[i][j][k] = 0;
373 HISTOGRAMS.hist_unlike_3d_coarse[i][j][k] = 0;
374 HISTOGRAMS.htmp_like_3d_fine[i][j][k] = 0;
375 HISTOGRAMS.htmp_unlike_3d_fine[i][j][k] = 0;
376 HISTOGRAMS.htmp_like_3d_coarse[i][j][k] = 0;
377 HISTOGRAMS.htmp_unlike_3d_coarse[i][j][k] = 0;
378 HISTOGRAMS.href_like_3d_fine[i][j][k] = 0;
379 HISTOGRAMS.href_unlike_3d_fine[i][j][k] = 0;
380 HISTOGRAMS.href_like_3d_coarse[i][j][k] = 0;
381 HISTOGRAMS.href_unlike_3d_coarse[i][j][k] = 0;
382 HISTOGRAMS.hinc_like_3d_fine[i][j][k] = 0;
383 HISTOGRAMS.hinc_unlike_3d_fine[i][j][k] = 0;
384 HISTOGRAMS.hinc_like_3d_coarse[i][j][k] = 0;
385 HISTOGRAMS.hinc_unlike_3d_coarse[i][j][k] = 0;
386 }
387/*********************************************/
388
389
390/*********************************************/
391
392// cout<<" FINISHED"<<endl;
393
394}
395
396
397/*****************************************************************************************/
398
399
400Int_t THBTprocessor::ImportParticles(TClonesArray *particles, Option_t *option)
401 {
402 //Copy particle data into TClonesArray
403 if (particles == 0) return 0;
404 TClonesArray &Particles = *particles;
405 Particles.Clear();
406
407 Int_t nrpart = 0;
408 for (Int_t i=0; i < TRK_MAXLEN; i++)
409 {
410
411
412 if (TRACK1.trk_E[i] == 0.) continue;
413
414 Float_t px = TRACK1.trk_px[i];
415 Float_t py = TRACK1.trk_py[i];
416 Float_t pz = TRACK1.trk_pz[i];
417// Float_t pE = TRACK.trk_E[i];
418 Float_t mass = PARTICLE.part_mass[TRACK1.trk_ge_pid[i]];
419
420 new(Particles[nrpart++]) TParticle(0,0,0,0,0,0,px,py,pz,
421 TMath::Sqrt(mass*mass+px*px+py*py+pz*pz),
422 0,0,0,0);
423 }
424 return nrpart;
425 }
426
427/*****************************************************************************************/
428
429void THBTprocessor::PrintEvent()
430 {
431 //Prints all particles (common block data)
432 cout<<"Print Event"<<endl;
433 for (Int_t i=0; i<TRK_MAXLEN;i++)
434 {
435 if(TRACK1.trk_E[i]==0.) continue;
436
437 cout<<"trk_id: "<<TRACK1.trk_id[i]<<" trk_px :"<<TRACK1.trk_px[i]<<" trk_py :"<<TRACK1.trk_py[i]<<" trk_pz :"<<TRACK1.trk_pz[i]<<endl;
438 cout<<" trk_E: "<<TRACK1.trk_E[i]<<" trk_pt: "<<TRACK1.trk_pt[i]<<" trk_phi: "<<TRACK1.trk_phi[i]<<" trk_eta: "<<TRACK1.trk_eta[i]<<endl;
439 }
440 }
441
442
443/*****************************************************************************************/
444void THBTprocessor::DumpSettings()
445{
446 //prints values set in common blocks
447 ctest();
448}