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