Coding conventions violations removed
[u/mrichter/AliRoot.git] / THbtp / THBTprocessor.cxx
1 #include "THBTprocessor.h"
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
22 #include <TParticle.h>
23 #include <TMath.h>
24
25 #include <Riostream.h>
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
37 ClassImp(THBTprocessor)
38  
39 extern "C" void  type_of_call hbtprocessor();   
40 extern "C" void  type_of_call ctest();
41
42
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
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  }
57 /*****************************************************************************************/
58
59 void THBTprocessor::GenerateEvent() const
60 {
61 //Starts processing
62
63   Info("GenerateEvent","Entering Fortran");
64   ctest();
65   hbtprocessor();
66   Info("GenerateEvent","Exited Fortran");
67   if(PARAMETERS.errorcode)
68     {
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());
75     }
76   else
77     Info("GenerateEvent","GOOD ! HBT Processor finished without errors");
78 }
79 /*****************************************************************************************/
80
81 void THBTprocessor::Initialize() const
82
83   //IT RESETS ALL THE PREVIOUS SETTINGS
84   //Call this method to set default values in PARAMETERS & MESH
85   //and zero other common block
86   
87   if(gDebug) 
88    Info("Initialize","Setting Default valuses in all COMMON BLOCKS");
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 /*********************************************/
185
186   Int_t i; //loop variable
187   
188   for (i =0; i<TRK_MAXLEN;i++)
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
213   for (i =0; i<TRK2_MAXLEN;i++)
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
238   for (i =0; i<SEC_MAXLEN;i++)
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
250   for (i =0; i<SEC_MAXLEN2;i++)
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
262   for (i =0; i<PART_MAXLEN;i++)
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 /*********************************************/
272   for (i =0; i<MAX_C2_1D;i++)
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 /*********************************************/
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++)
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
328  for (i =0; i<MAX_EVENTS;i++) 
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 /*********************************************/
341  for (i =0; i<MAX_C2_COUL;i++) 
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 /*********************************************/
349  for (i =0; i<MAX_H_1D;i++) 
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
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++) 
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
419 Int_t THBTprocessor::ImportParticles(TClonesArray *particles, Option_t */*option*/)
420  {
421   //Copy particle data into TClonesArray
422   if (particles == 0) return 0;
423   TClonesArray &rparticles = *particles;
424   rparticles.Clear();
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     
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),
441                                          0,0,0,0);
442    }
443   return nrpart;        
444  }
445
446 /*****************************************************************************************/
447
448 void THBTprocessor::PrintEvent() const
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 /*****************************************************************************************/
463 void THBTprocessor::DumpSettings() const
464 {
465  //prints values set in common blocks
466   ctest();
467 }