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 }