THbtp/THBTprocessor.cxx

   1 #include <TGenerator.h>
   2 #include "THBTprocessor.h"
   3 #include <TParticle.h>
   4 #include <TMath.h>
   5
   6 #include <iostream.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
  18 ClassImp(THBTprocessor)
  19
  20 extern "C" void  type_of_call hbtprocessor();
  21 extern "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 /*************************************************/
  30 THBTprocessor::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
  44 void 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 /*****************************************************************************************/
  63 void 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
 400 Int_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
 429 void 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 /*****************************************************************************************/