9d9cd6117e4b157cbafb657c7c9cc5a19b8bf896
[u/mrichter/AliRoot.git] / THbtp / THBTprocessor.cxx
1
2 #include <TGenerator.h>
3 #include "THBTprocessor.h"
4 #include <TParticle.h>
5 #include <TMath.h>
6
7 #include <iostream.h>
8 #ifndef WIN32
9 # define hbtprocessor hbtprocessor_
10 # define ctest ctest_
11 # define type_of_call
12 #else
13 # define hbtprocessor HBTPROCESSOR
14 # define ctest CTEST
15 # define type_of_call _stdcall
16 #endif    
17
18
19 ClassImp(THBTprocessor)
20  
21 extern "C" void  type_of_call hbtprocessor();   
22 extern "C" void  type_of_call ctest();
23
24 //Wrapper class to HBT processor fortran code.
25 //For more information see AliGenHBTprocessor class
26 //HBT processor is written by Lanny Ray
27 //Wrapper class written by Piotr Krzysztof Skowronski (Piotr.Skowronski@cern.ch)
28 //
29
30 /*************************************************/
31 THBTprocessor::THBTprocessor()// it is better not to intialize it:TGenerator("THBTprocessor","THBTprocessor")
32                               //because it allocates memmory for TObjArray::fParticles which is not used in our case
33                               //and we are using TClonesArray in import paerticles
34  {
35  //constructor
36   PARAMETERS.ALICE            = 1; //flag that we are working in AliRoot (0==STAR stand alone mode)
37   PARAMETERS.pi               = TMath::Pi();//3.141592654;
38   PARAMETERS.rad              = 180.0/TMath::Pi();
39   PARAMETERS.hbc              = 0.19732891;
40   
41   fParticles = 0;
42   Initialize(); //Enforce initialization (cleaning all commons)
43  }
44
45 void THBTprocessor::GenerateEvent()
46 {
47 //Starts processing
48
49   cout<<endl<<"Entering Fortran"<<endl;
50   ctest();
51   hbtprocessor();
52   cout<<endl<<"Exited Fortran"<<endl;
53   if(PARAMETERS.errorcode)
54     {
55       Fatal("THBTprocessor::GenerateEvent()","HBT Processor (fortran part) finished with errors\n \
56              See hbt_simulation.out file for more detailed information");
57     }
58   else
59       cout<<endl<<"GOOD ! HBT Processor finished without errors"<<endl;
60   cout<<"PARAMETERS.errorcode= "<<PARAMETERS.errorcode<<endl;
61   
62 }
63 /*****************************************************************************************/
64 void THBTprocessor::Initialize()
65 {
66   //Call this method to set default values in PARAMETERS & MESH
67   //and zero other common block
68   
69   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   for (Int_t i =0; i<TRK_MAXLEN;i++)
168    {
169     TRACK1.trk_id[i] = 0;
170     TRACK1.trk_px_sec[i] = 0; 
171     TRACK1.trk_py_sec[i] = 0;
172     TRACK1.trk_pz_sec[i] = 0;
173     TRACK1.trk_sector[i] = 0;
174     TRACK1.trk_flag[i] = 0;
175     TRACK1.trk_out_flag[i] = 0;
176     TRACK1.trk_merge_flag[i] = 0;
177     TRACK1.trk_ge_pid[i] = 0;
178     TRACK1.trk_start_vertex[i] = 0;
179     TRACK1.trk_stop_vertex[i] = 0;
180     TRACK1.trk_event_line[i] = 0;
181     TRACK1.trk_px[i] = 0;
182     TRACK1.trk_py[i] = 0;
183     TRACK1.trk_pz[i] = 0;
184     TRACK1.trk_E[i] = 0;
185     TRACK1.trk_pt[i] = 0;
186     TRACK1.trk_phi[i] = 0;
187     TRACK1.trk_eta[i] = 0; 
188    }
189   
190 /*********************************************/
191
192   for (Int_t i =0; i<TRK2_MAXLEN;i++)
193    {
194     TRACK2.trk2_id[i] = 0;
195     TRACK2.trk2_px_sec[i] = 0; 
196     TRACK2.trk2_py_sec[i] = 0;
197     TRACK2.trk2_pz_sec[i] = 0;
198     TRACK2.trk2_sector[i] = 0;
199     TRACK2.trk2_flag[i] = 0;
200     TRACK2.trk2_out_flag[i] = 0;
201     TRACK2.trk2_merge_flag[i] = 0;
202     TRACK2.trk2_ge_pid[i] = 0;
203     TRACK2.trk2_start_vertex[i] = 0;
204     TRACK2.trk2_stop_vertex[i] = 0;
205     TRACK2.trk2_event_line[i] = 0;
206     TRACK2.trk2_px[i] = 0;
207     TRACK2.trk2_py[i] = 0;
208     TRACK2.trk2_pz[i] = 0;
209     TRACK2.trk2_E[i] = 0;
210     TRACK2.trk2_pt[i] = 0;
211     TRACK2.trk2_phi[i] = 0;
212     TRACK2.trk2_eta[i] = 0; 
213    }
214
215 /*********************************************/
216
217   for (Int_t i =0; i<SEC_MAXLEN;i++)
218    {
219     SEC_TRK_MAP.stm_sec_id [i] = 0;
220     SEC_TRK_MAP.stm_n_trk_sec[i] = 0; 
221     SEC_TRK_MAP.stm_flag[i] = 0;
222     
223     for (Int_t j=0; j<MAX_TRK_SEC;j++)
224        SEC_TRK_MAP.stm_track_id[i][j] = 0;
225    }
226
227 /*********************************************/
228
229   for (Int_t i =0; i<SEC_MAXLEN2;i++)
230    {
231     SEC_TRK_MAP2.stm_sec_id2[i] = 0;
232     SEC_TRK_MAP2.stm_n_trk_sec2[i] = 0; 
233     SEC_TRK_MAP2.stm_flag2[i] = 0;
234     
235     for (Int_t j=0; j<MAX_TRK_SEC;j++)
236        SEC_TRK_MAP2.stm_track_id2[i][j] = 0;
237    }
238
239 /*********************************************/
240
241   for (Int_t i =0; i<PART_MAXLEN;i++)
242    {
243      PARTICLE.part_id[i] = 0;
244      PARTICLE.part_charge[i] = 0;
245      PARTICLE.part_mass[i] = 0;
246      PARTICLE.part_lifetime[i] = 0;
247    }
248
249
250 /*********************************************/
251   for (Int_t i =0; i<MAX_C2_1D;i++)
252    {
253      CORRELATIONS.c2mod_like_1d[i] = 0; 
254      CORRELATIONS.c2mod_unlike_1d[i] = 0;
255      CORRELATIONS.c2fit_like_1d[i] = 0;
256      CORRELATIONS.c2fit_unlike_1d[i] = 0;
257      CORRELATIONS.c2err_like_1d[i] = 0;
258      CORRELATIONS.c2err_unlike_1d[i] = 0;
259    }
260 /*********************************************/
261   for (Int_t i =0; i<MAX_C2_3D;i++)
262    for (Int_t j =0; j<MAX_C2_3D;j++)
263     for (Int_t k =0; k<MAX_C2_3D;k++)
264      {
265       CORRELATIONS.c2mod_like_3d_fine[i][j][k] = 0.0;
266       CORRELATIONS.c2mod_unlike_3d_fine[i][j][k] = 0.0;
267       CORRELATIONS.c2mod_like_3d_coarse[i][j][k] = 0.0;
268       CORRELATIONS.c2mod_unlike_3d_coarse[i][j][k] = 0.0;
269       CORRELATIONS.c2fit_like_3d_fine[i][j][k] = 0.0;
270       CORRELATIONS.c2fit_unlike_3d_fine[i][j][k] = 0.0;
271       CORRELATIONS.c2fit_like_3d_coarse[i][j][k] = 0.0;
272       CORRELATIONS.c2fit_unlike_3d_coarse[i][j][k] = 0.0;
273       CORRELATIONS.c2err_like_3d_fine[i][j][k] = 0.0;
274       CORRELATIONS.c2err_unlike_3d_fine[i][j][k] = 0.0;
275       CORRELATIONS.c2err_like_3d_coarse[i][j][k] = 0.0;
276       CORRELATIONS.c2err_unlike_3d_coarse[i][j][k] = 0.0;
277      }
278 /*********************************************/
279
280    EVENT_SUMMARY.niter_mean = 0.0; 
281    EVENT_SUMMARY.niter_rms = 0.0;
282  
283    EVENT_SUMMARY.npart1_mean = 0.0; 
284    EVENT_SUMMARY.npart1_rms = 0.0;
285  
286    EVENT_SUMMARY.npart2_mean = 0.0; 
287    EVENT_SUMMARY.npart2_rms = 0.0;
288  
289    EVENT_SUMMARY.npart_tot_mean = 0.0; 
290    EVENT_SUMMARY.npart_tot_rms = 0.0;
291  
292    EVENT_SUMMARY.nsec_flag_mean = 0.0; 
293    EVENT_SUMMARY.nsec_flag_rms = 0.0;
294  
295    EVENT_SUMMARY.frac_trks_out_mean = 0.0;
296    EVENT_SUMMARY.frac_trks_out_rms = 0.0;
297  
298    EVENT_SUMMARY.frac_trks_flag_mean = 0.0;
299    EVENT_SUMMARY.frac_trks_flag_rms = 0.0;
300  
301    EVENT_SUMMARY.chi_l1d_mean = 0.0;
302    EVENT_SUMMARY.chi_l1d_rms = 0.0;
303
304    EVENT_SUMMARY.chi_u1d_mean = 0.0;
305    EVENT_SUMMARY.chi_u1d_rms = 0.0;
306
307  for (Int_t i =0; i<MAX_EVENTS;i++) 
308   {
309     EVENT_SUMMARY.n_part_used_1_store[i] = 0.0; 
310     EVENT_SUMMARY.n_part_used_2_store[i] = 0.0; 
311     EVENT_SUMMARY.n_part_tot_store[i] = 0.0; 
312     EVENT_SUMMARY.num_sec_flagged_store[i] = 0.0; 
313     EVENT_SUMMARY.frac_trks_out[i] = 0.0; 
314     EVENT_SUMMARY.frac_trks_flag[i] = 0.0;
315     EVENT_SUMMARY.chisq_like_1d_store[i] = 0.0;         
316     EVENT_SUMMARY.num_iter[i] = 0.0; 
317     EVENT_SUMMARY.chisq_unlike_1d_store[i] = 0.0;       
318   } 
319 /*********************************************/
320  for (Int_t i =0; i<MAX_C2_COUL;i++) 
321   {
322       COULMB.c2_coul_like[i] = 0.0;
323       COULMB.c2_coul_unlike[i] = 0.0;
324       COULMB.q_coul[i] = 0.0; 
325  
326   }
327 /*********************************************/
328  for (Int_t i =0; i<MAX_H_1D;i++) 
329    {  
330       HISTOGRAMS.hist1_pt_1[i] = 0;
331       HISTOGRAMS.hist1_phi_1[i] = 0;
332       HISTOGRAMS.hist1_eta_1[i] = 0;
333       HISTOGRAMS.hist1_pt_2[i] = 0;
334       HISTOGRAMS.hist1_phi_2[i] = 0;
335       HISTOGRAMS.hist1_eta_2[i] = 0;
336       HISTOGRAMS.htmp1_pt_1[i] = 0;
337       HISTOGRAMS.htmp1_phi_1[i] = 0;
338       HISTOGRAMS.htmp1_eta_1[i] = 0;
339       HISTOGRAMS.htmp1_pt_2[i] = 0;
340       HISTOGRAMS.htmp1_phi_2[i] = 0;
341       HISTOGRAMS.htmp1_eta_2[i] = 0;
342       HISTOGRAMS.href1_pt_1[i] = 0;
343       HISTOGRAMS.href1_phi_1[i] = 0;
344       HISTOGRAMS.href1_eta_1[i] = 0;
345       HISTOGRAMS.href1_pt_2[i] = 0;
346       HISTOGRAMS.href1_phi_2[i] = 0;
347       HISTOGRAMS.href1_eta_2[i] = 0;
348       HISTOGRAMS.hinc1_pt_1[i] = 0;
349       HISTOGRAMS.hinc1_phi_1[i] = 0;
350       HISTOGRAMS.hinc1_eta_1[i] = 0;
351       HISTOGRAMS.hinc1_pt_2[i] = 0;
352       HISTOGRAMS.hinc1_phi_2[i] = 0;
353       HISTOGRAMS.hinc1_eta_2[i] = 0;
354       HISTOGRAMS.hist_like_1d[i] = 0;
355       HISTOGRAMS.hist_unlike_1d[i] = 0;
356       HISTOGRAMS.htmp_like_1d[i] = 0;
357       HISTOGRAMS.htmp_unlike_1d[i] = 0;
358       HISTOGRAMS.href_like_1d[i] = 0;
359       HISTOGRAMS.href_unlike_1d[i] = 0;
360       HISTOGRAMS.hinc_like_1d[i] = 0;
361       HISTOGRAMS.hinc_unlike_1d[i] = 0;
362   }
363
364  for (Int_t i =0; i<MAX_H_3D;i++) 
365    for (Int_t j =0; j<MAX_H_3D;j++) 
366      for (Int_t k =0; k<MAX_H_3D;k++) 
367        {
368           HISTOGRAMS.hist_like_3d_fine[i][j][k] = 0;
369           HISTOGRAMS.hist_unlike_3d_fine[i][j][k] = 0;
370           HISTOGRAMS.hist_like_3d_coarse[i][j][k] = 0;
371           HISTOGRAMS.hist_unlike_3d_coarse[i][j][k] = 0;
372           HISTOGRAMS.htmp_like_3d_fine[i][j][k] = 0;
373           HISTOGRAMS.htmp_unlike_3d_fine[i][j][k] = 0;
374           HISTOGRAMS.htmp_like_3d_coarse[i][j][k] = 0;
375           HISTOGRAMS.htmp_unlike_3d_coarse[i][j][k] = 0;
376           HISTOGRAMS.href_like_3d_fine[i][j][k] = 0;
377           HISTOGRAMS.href_unlike_3d_fine[i][j][k] = 0;
378           HISTOGRAMS.href_like_3d_coarse[i][j][k] = 0;
379           HISTOGRAMS.href_unlike_3d_coarse[i][j][k] = 0;
380           HISTOGRAMS.hinc_like_3d_fine[i][j][k] = 0;
381           HISTOGRAMS.hinc_unlike_3d_fine[i][j][k] = 0;
382           HISTOGRAMS.hinc_like_3d_coarse[i][j][k] = 0;
383           HISTOGRAMS.hinc_unlike_3d_coarse[i][j][k] = 0;
384        }
385 /*********************************************/
386
387
388 /*********************************************/
389
390 //  cout<<" FINISHED"<<endl;
391   
392 }
393
394
395 /*****************************************************************************************/
396
397
398 Int_t THBTprocessor::ImportParticles(TClonesArray *particles, Option_t *option)
399  {
400   //Copy particle data into TClonesArray
401   if (particles == 0) return 0;
402   TClonesArray &Particles = *particles;
403   Particles.Clear();
404  
405   Int_t nrpart = 0;
406   for (Int_t i=0; i < TRK_MAXLEN; i++) 
407    {
408    
409     
410       if (TRACK1.trk_E[i] == 0.) continue;
411     
412       Float_t px   = TRACK1.trk_px[i];
413       Float_t py   = TRACK1.trk_py[i];
414       Float_t pz   = TRACK1.trk_pz[i];
415 //    Float_t pE   = TRACK.trk_E[i];
416       Float_t mass = PARTICLE.part_mass[TRACK1.trk_ge_pid[i]];
417     
418       new(Particles[nrpart++]) TParticle(0,0,0,0,0,0,px,py,pz,
419                                          TMath::Sqrt(mass*mass+px*px+py*py+pz*pz),
420                                          0,0,0,0);
421    }
422   return nrpart;        
423  }
424
425 /*****************************************************************************************/
426
427 void THBTprocessor::PrintEvent()
428  {
429  //Prints all particles (common block data)  
430    cout<<"Print Event"<<endl;
431    for (Int_t i=0; i<TRK_MAXLEN;i++)
432     {
433       if(TRACK1.trk_E[i]==0.) continue;
434
435       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;
436       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;
437     }
438  }
439
440
441 /*****************************************************************************************/