1 #include "THBTprocessor.h"
2 //_____________________________________________________________________________
3 ///////////////////////////////////////////////////////////////////////////////
5 // class THBTprocessor //
7 // Wrapper class to HBT processor fortran code. //
8 // For more information see AliGenHBTprocessor class //
9 // HBT processor is written by Lanny Ray //
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 //
16 // Wrapper class written by //
17 // Piotr Krzysztof Skowronski (Piotr.Skowronski@cern.ch) //
19 ///////////////////////////////////////////////////////////////////////////////
22 #include <TParticle.h>
25 #include <Riostream.h>
27 # define hbtprocessor hbtprocessor_
31 # define hbtprocessor HBTPROCESSOR
33 # define type_of_call _stdcall
37 ClassImp(THBTprocessor)
39 extern "C" void type_of_call hbtprocessor();
40 extern "C" void type_of_call ctest();
43 /*************************************************/
44 THBTprocessor::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
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;
55 Initialize(); //Enforce initialization (cleaning all commons)
57 /*****************************************************************************************/
59 void THBTprocessor::GenerateEvent() const
63 Info("GenerateEvent","Entering Fortran");
66 Info("GenerateEvent","Exited Fortran");
67 if(PARAMETERS.errorcode)
69 TString message("HBT Processor (fortran part) finished with errors\n");
70 message+="Error code is ";
71 message+=PARAMETERS.errorcode;
73 message+="See hbt_simulation.out file for more detailed information.";
74 Fatal("GenerateEvent","%s",message.Data());
77 Info("GenerateEvent","GOOD ! HBT Processor finished without errors");
79 /*****************************************************************************************/
81 void THBTprocessor::Initialize() const
83 //IT RESETS ALL THE PREVIOUS SETTINGS
84 //Call this method to set default values in PARAMETERS & MESH
85 //and zero other common block
88 Info("Initialize","Setting Default valuses in all COMMON BLOCKS");
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;
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;
112 PARAMETERS.Rout = 7.0;
113 PARAMETERS.Rlong = 4.0;
114 PARAMETERS.Rperp = 6.0;
115 PARAMETERS.Rparallel = 4.0;
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;
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;
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;
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;
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
153 /*********************************************/
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
184 /*********************************************/
186 Int_t i; //loop variable
188 for (i =0; i<TRK_MAXLEN;i++)
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;
206 TRACK1.trk_pt[i] = 0;
207 TRACK1.trk_phi[i] = 0;
208 TRACK1.trk_eta[i] = 0;
211 /*********************************************/
213 for (i =0; i<TRK2_MAXLEN;i++)
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;
236 /*********************************************/
238 for (i =0; i<SEC_MAXLEN;i++)
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;
244 for (Int_t j=0; j<MAX_TRK_SEC;j++)
245 SEC_TRK_MAP.stm_track_id[i][j] = 0;
248 /*********************************************/
250 for (i =0; i<SEC_MAXLEN2;i++)
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;
256 for (Int_t j=0; j<MAX_TRK_SEC;j++)
257 SEC_TRK_MAP2.stm_track_id2[i][j] = 0;
260 /*********************************************/
262 for (i =0; i<PART_MAXLEN;i++)
264 PARTICLE.part_id[i] = 0;
265 PARTICLE.part_charge[i] = 0;
266 PARTICLE.part_mass[i] = 0;
267 PARTICLE.part_lifetime[i] = 0;
271 /*********************************************/
272 for (i =0; i<MAX_C2_1D;i++)
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;
281 /*********************************************/
282 for (i =0; i<MAX_C2_3D;i++)
283 for (Int_t j =0; j<MAX_C2_3D;j++)
284 for (Int_t k =0; k<MAX_C2_3D;k++)
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;
299 /*********************************************/
301 EVENT_SUMMARY.niter_mean = 0.0;
302 EVENT_SUMMARY.niter_rms = 0.0;
304 EVENT_SUMMARY.npart1_mean = 0.0;
305 EVENT_SUMMARY.npart1_rms = 0.0;
307 EVENT_SUMMARY.npart2_mean = 0.0;
308 EVENT_SUMMARY.npart2_rms = 0.0;
310 EVENT_SUMMARY.npart_tot_mean = 0.0;
311 EVENT_SUMMARY.npart_tot_rms = 0.0;
313 EVENT_SUMMARY.nsec_flag_mean = 0.0;
314 EVENT_SUMMARY.nsec_flag_rms = 0.0;
316 EVENT_SUMMARY.frac_trks_out_mean = 0.0;
317 EVENT_SUMMARY.frac_trks_out_rms = 0.0;
319 EVENT_SUMMARY.frac_trks_flag_mean = 0.0;
320 EVENT_SUMMARY.frac_trks_flag_rms = 0.0;
322 EVENT_SUMMARY.chi_l1d_mean = 0.0;
323 EVENT_SUMMARY.chi_l1d_rms = 0.0;
325 EVENT_SUMMARY.chi_u1d_mean = 0.0;
326 EVENT_SUMMARY.chi_u1d_rms = 0.0;
328 for (i =0; i<MAX_EVENTS;i++)
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;
340 /*********************************************/
341 for (i =0; i<MAX_C2_COUL;i++)
343 COULMB.c2_coul_like[i] = 0.0;
344 COULMB.c2_coul_unlike[i] = 0.0;
345 COULMB.q_coul[i] = 0.0;
348 /*********************************************/
349 for (i =0; i<MAX_H_1D;i++)
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;
385 for (i =0; i<MAX_H_3D;i++)
386 for (Int_t j =0; j<MAX_H_3D;j++)
387 for (Int_t k =0; k<MAX_H_3D;k++)
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;
406 /*********************************************/
409 /*********************************************/
411 // cout<<" FINISHED"<<endl;
416 /*****************************************************************************************/
419 Int_t THBTprocessor::ImportParticles(TClonesArray *particles, Option_t */*option*/)
421 //Copy particle data into TClonesArray
422 if (particles == 0) return 0;
423 TClonesArray &rparticles = *particles;
427 for (Int_t i=0; i < TRK_MAXLEN; i++)
431 if (TRACK1.trk_E[i] == 0.) continue;
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]];
439 new(rparticles[nrpart++]) TParticle(0,0,0,0,0,0,px,py,pz,
440 TMath::Sqrt(mass*mass+px*px+py*py+pz*pz),
446 /*****************************************************************************************/
448 void THBTprocessor::PrintEvent() const
450 //Prints all particles (common block data)
451 cout<<"Print Event"<<endl;
452 for (Int_t i=0; i<TRK_MAXLEN;i++)
454 if(TRACK1.trk_E[i]==0.) continue;
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;
462 /*****************************************************************************************/
463 void THBTprocessor::DumpSettings() const
465 //prints values set in common blocks