HBTP/hbt_event_processor.f

   1 CCC ********************************************************************
   2 CCC  Modifications for ALICE-ROOT application at CERN - 12/15/2000
   3 CCC  1.  Removed all Fortran Data Structures in favor of labelled
   4 CCC      common blocks.  The syntax of the structure_variable to
   5 CCC      common variable name change is:
   6 CCC      In STAR code    In ALICE code
   7 CCC        A(i).B     ==>  A_B(i)
   8 CCC        A(i).B(j)  ==>  A_B(j,i)
   9 CCC  2.  All remaining references in the comments and write statements
  10 CCC      to the data structures are interpreted as applying to the
  11 CCC      new common variables.
  12 CCC  3.  The UNIX random number generator, ran(), was replaced with
  13 CCC      a function which calls a modified version of the CERNLIB
  14 CCC      random number generator, ranlux, herein called ranlux2.  The
  15 CCC      latter allows a user supplied seed value.
  16 CCC  4.  Increased the following array sizes for the LHC Pb+Pb design
  17 CCC      multiplicity criteria of dN_{ch}/dy = 8000 assuming 80% of
  18 CCC      this is pi(+) and pi(-), which are the largest populations
  19 CCC      that this code is required to process.
  20 CCC
  21 CCC      The following are for the max number of tracks that can be
  22 CCC      stored in each sector without overflow:
  23 CCC         common_mesh.inc - increased max_trk_save from 30 to 150
  24 CCC         common_sec_track.inc  - inc. max_trk_sec from 30 to 150
  25 CCC         common_sec_track2.inc - inc. max_trk_sec2 from 30 to 150
  26 CCC
  27 CCC      The following determine the maximum number of tracks that can
  28 CCC      be processed in an event:
  29 CCC         common_track.inc - increased trk_maxlen from 6500 to 25000
  30 CCC         common_track2.inc - increased trk2_maxlen from 6500 to 25000
  31 CCC
  32 C
  33 C     DESCRIPTION OF METHOD:
  34 C     =====================
  35 C
  36 C          This program produces relativistic heavy-ion collision
  37 C     events in which the particle momenta for selected particle ID
  38 C     types and for selected kinematic acceptance ranges are randomly
  39 C     adjusted such that specified one-body distributions and two-body
  40 C     correlation functions are reproduced.  The input to the code may
  41 C     be a set of events from any STAR event generator, so long as the
  42 C     format is in the STAR Geant (GSTAR) text format standard (see
  43 C     STAR Note #235).  The basic method is similar to that of Metropolis
  44 C     et al. and is fully described in Ray and Hoffmann, Phys. Rev. C 54,
  45 C     2582 (1996).  Briefly the steps in the algorithm include:
  46 C
  47 C        (1) For an initial particle distribution of specified particle
  48 C            ID types (maximum of two types allowed) and momentum
  49 C            acceptance ranges [given in terms of transverse momentum
  50 C            (p_T), azimuthal angle (phi) and pseudorapidity (eta)] the
  51 C            momentum vector of one particle is randomly shifted within
  52 C            a specified range from its initial value.  The shifts are
  53 C            done for px, py and pz independently.
  54 C
  55 C        (2) New one-body and two-body histograms, as well as the two-body
  56 C            correlation function are calculated.
  57 C
  58 C        (3) If the random momentum shift results in an improved overall
  59 C            chi-square, obtained by comparison with a specified reference
  60 C            for the one-body distribution and the two-body correlation
  61 C            model, then the new momentum vector is retained.  If not,
  62 C            then the vector is restored to its starting value.
  63 C
  64 C        (4) Steps 1-3 are repeated for each accepted particle in the
  65 C            event.
  66 C
  67 C        (5) The entire process, steps 1-4, is repeated until either a
  68 C            satisfactory fit to the model distributions are obtained or
  69 C            the maximum number of iterations is reached.
  70 C
  71 C        (6) Once the iterative process is complete, the input event file
  72 C            is copied directly to an output event file where the adjusted
  73 C            momentum values for the accepted tracks replace that in the
  74 C            input file.  The event output file is in the GSTAR standard
  75 C            text format.  This event output file may be processed again
  76 C            by this code in order to generate correlations for other
  77 C            particle types or for different kinematic ranges.  The file
  78 C            is suitable for input into the STAR version of Geant-3, called
  79 C            GSTAR (STAR Note 235).
  80 C
  81 C     In order to reduce cpu demand the particle momenta are sorted into
  82 C     momentum space cells or sectors.  In forming particle pairs only those
  83 C     particles in the same, or adjacent cells are used.  For large events
  84 C     this vastly reduces the required cpu time.  The penalty is that the
  85 C     coding becomes more complicated.  Much of the present code is devoted
  86 C     to the necessary bookeeping chores that must be done in order to
  87 C     determine which cell the tracks are in and if they move to new cells
  88 C     during the track adjustment procedure.  Information about the
  89 C     momentum space cells are contained in the data structure /sec_trk_map/.
  90 C
  91 C     The sector size must therefore be scaled with the specified correlation
  92 C     range.  All particles will be paired with all possible partners out
  93 C     to Q's equal to the smallest dimension of the momentum space sectors.
  94 C     Particle pairs with Q greater than this sector dimension will suffer
  95 C     reduced acceptance, finally being completely cut-off for Q ~> 2 times
  96 C     the diagonal length thru a sector.
  97 C
  98 C          In order to generate momentum correlations for particle types
  99 C     having low multiplicity it is necessary for the user to supply this
 100 C     code with an artificially enhanced multiplicity along with a track-
 101 C     write-output fractional acceptance factor (see input variable
 102 C     'trk_accep').  For example, if the user wants to generate HBT
 103 C     correlations for K0-shorts but the assumed multiplicity is too
 104 C     low for the present algorithm to work, the user may increase the
 105 C     input K0-short multiplicity, for instance by a factor of 5, then
 106 C     run the code and set trk_accep = 1/5 = 0.2 in order to randomly
 107 C     reject about 80% of the K0-shorts.  The track rejection is done
 108 C     after the track adjustment procedure is completed.  This procedure
 109 C     should preserve the built-in correlations in the final output
 110 C     event file, although with reduced statistics.
 111 C          Another approach for handling low multiplicity events would
 112 C     be to combine particles from several events, carry out the track
 113 C     adjustment procedure, then separate the tracks into their original
 114 C     events.  This method must insure that no bias is included due to
 115 C     the order of processing the tracks from the first, second, etc.
 116 C     events.  This latter method, once proven, could be used for
 117 C     the low multiplicity particles from any event generator.  For
 118 C     the present version of the code the low multiplicity HBT studies
 119 C     must utilize a Monte Carlo multiplicity generator.
 120 C
 121 C          The code may also be used to read in a previously correlated
 122 C     set of events and either compute the reference histograms or read in
 123 C     the references, and then compute the correlations for each event and
 124 C     the inclusive correlations without doing the track momentum adjustment
 125 C     procedure.  This feature may be used, for example, to study the
 126 C     correlations that result in one set of coordinates for events fitted
 127 C     to correlations with respect to a different set of coordinates.  For
 128 C     example, fit the correlations to the Y-K-P form and then evaluate
 129 C     the side-out-long correlations, or vice-versa.
 130 C
 131 C     TWO-BODY REFERENCE HISTOGRAMS:
 132 C     =============================
 133 C
 134 C          In order to calculate the correlations, an uncorrelated two-body
 135 C     reference spectrum is needed.  The program will calculate this
 136 C     quantity by forming pairs of particles from different events in the
 137 C     input file.  For the particle ID type(s) and momentum acceptance
 138 C     the code forms all possible pairs (given the cell substructure) by
 139 C     mixing particles from event#1 with those in event#2, then particles
 140 C     from event#2 are mixed with particles from event#3, then events 3
 141 C     and 4 are mixed, and so on.  In this way ample statistics may be
 142 C     achieved for the reference distributions.  These reference distributions
 143 C     can be written out to file and re-used in subsequent runs.  Since
 144 C     all events in the input event file are used in generating the
 145 C     reference distribution, it is imperative that these events be physically
 146 C     similar.
 147 C
 148 C     ONE-BODY REFERENCE HISTOGRAMS:
 149 C     =============================
 150 C
 151 C          Inclusive sums of the accepted particles for all events in the
 152 C     input event file determine the one-body reference distributions.
 153 C     These are used to constrain the momentum vector shifts.  Although
 154 C     the one-body distributions from realistic event generators are fully
 155 C     three-dimensional, the present code is restricted to only work with
 156 C     the one-dimensional projections onto p_T, phi and eta.  In other words,
 157 C     the p_T distribution used in this code is formed by integrating
 158 C     the particle distributions over (phi,eta) over the momentum acceptance
 159 C     range.  No particle distribution models are built into the code.
 160 C     The one-body reference distributions are either read-in or determined
 161 C     from the events in the input event text file.
 162 C
 163 C     TWO-BODY CORRELATION MODELS:
 164 C     ===========================
 165 C
 166 C          The code permits both 1-dimensional and 3-dimensional two-body
 167 C     correlation models.  These may be fitted separately or simultaneously.
 168 C     The source may include a mixture of incoherent and coherent components;
 169 C     Coulomb corrections are also included.  The general form assumed
 170 C     [see Cramer and Kadija, Phys. Rev. C 53, 908 (1996)] is:
 171 C
 172 C       C2 = 1 + lambda*(b**2) + 2.0*sqrt(lambda)*[1 - sqrt(lambda)]*b
 173 C
 174 C     where lambda is the usual chaoticity parameter.  The third term in
 175 C     this equation may be turned on or off.  Values of lambda < 1.0 may
 176 C     be used without the third term being included.  For 1-dimensional
 177 C     functions b is given by:
 178 C
 179 C       b = exp(- Q**2 * R**2 / 2)
 180 C
 181 C     where Q is either the invariant 4-momentum difference, the total
 182 C     4-momentum difference (i.e. time-like + space-like) or the
 183 C     3-vector momentum difference.  The 3-dimensional functions may be
 184 C     of the Bertsch-Pratt ``side-out-longitudinal'' form given by:
 185 C
 186 C       b = exp[(- Qside**2 * Rside**2 - Qout**2 * Rout**2
 187 C                - Qlong**2 * Rlong**2)/2]
 188 C
 189 C     where the ``out-long'' cross term is omitted.  The 3D function may
 190 C     also be in the Yano-Koonin-Podgoretski (YKP) form given by (for
 191 C     pairs in the A+A c.m. frame):
 192 C
 193 C       b = exp[(- Qperp**2 * Rperp**2 - Qparallel**2 * Rparallel**2
 194 C                - Qtime**2 * Rtime**2)/2]
 195 C
 196 C     where
 197 C                Qperp = transverse momentum difference
 198 C                Qparallel = Qlong = p_{1z} - p_{2z}
 199 C                Qtime = E_1 - E_2
 200 C
 201 C     The Coulomb correction may be omitted, or included in one of 3 ways:
 202 C
 203 C          (1) Point source Gamow factor
 204 C          (2) Finite source NA35 model (see Eq.(5) in Z. Phys. C73, 443
 205 C              (1997)) where
 206 C
 207 C              Coulomb correction = 1 + [G(q) - 1]*exp(-q/Q0)
 208 C
 209 C              and G(q) is the Gamow factor and q is the 3-momentum
 210 C              vector difference.
 211 C          (3) Finite source, Pratt integrated Coulomb wave function
 212 C              integrals, interpolated for a specific source radius
 213 C              from John Cramer's tables for discrete source radii.
 214 C
 215 C     These Coulomb correction factors multiply the above correlation
 216 C     function to give the total correlation function to be fitted for
 217 C     charged, like pairs. For charged, unlike pairs only the Coulomb
 218 C     (attractive) correlation function is used.
 219 C
 220 C     BINNING:
 221 C     =======
 222 C
 223 C          Several types of binning are done in the code. The one-body
 224 C     distributions are binned in p_t, phi and eta.  The full momentum
 225 C     space is subdivided into cells or sectors.  The 1D and 3D two-body
 226 C     distributions are binned into fine and coarse bins, where the fine
 227 C     bins occur at the smaller momentum difference range and the coarse
 228 C     bins at the larger.  For the 3D distributions the (1,1,1) coarse
 229 C     bin must coincide with the 3D fine bins.
 230 C
 231 C     SUMMARY OF EXTERNAL FILES:
 232 C     =========================
 233 C
 234 C File Unit#    File Name                Description
 235 C ----------------------------------------------------------------------------
 236 C     1      hbt_parameters.in          Input switches and controls
 237 C     2      event_text.in              Event generator output in GSTAR text
 238 C     3      event_line.flags           Generated tmp file, input line flags
 239 C     4      event_tracks.select        Generated tmp file, accep. tracks flg.
 240 C     7      hbt_log.out                Generated log file - error reports
 241 C     8      hbt_simulation.out         Generated main output file
 242 C     9      hbt_pair_reference.hist    Generated pair ref. histograms
 243 C     10     event_hbt_text.out         Gen. correlated event text file
 244 C     11     hbt_singles_reference.hist Gen. one-body ref. histograms
 245 C     12     event_text_aux.in          Tmp copy of event_text.in per event
 246 C     14     event_tracks_aux.select    Tmp copy of event_tracks.select/event
 247 C     21-27  cpp_*.dat (*=06,08...18)   Like pair Pratt Coulomb corrections.
 248 C     31-37  cpm_*.dat (*=06,08...18)   Unlike pair Pratt Coulomb corrects.
 249 C ----------------------------------------------------------------------------
 250 C
 251 C     Source of Data for ALICE Application:
 252 C     ====================================
 253 C
 254 C File Unit#    File Name                For ALICE File or Struc?
 255 C ----------------------------------------------------------------------------
 256 C     1      hbt_parameters.in          Call AliHbtp_ function
 257 C     2      event_text.in              Call AliHbtp_ function
 258 C     3      event_line.flags           File not used
 259 C     4      event_tracks.select        File not used
 260 C     7      hbt_log.out                File used as is
 261 C     8      hbt_simulation.out         File used as is
 262 C     9      hbt_pair_reference.hist    File used as is
 263 C     10     event_hbt_text.out         Call AliHbtp_ function
 264 C     11     hbt_singles_reference.hist File used as is
 265 C     12     event_text_aux.in          File not used
 266 C     14     event_tracks_aux.select    File not used
 267 C     21-27  cpp_*.dat (*=06,08...18)   Files are used as is
 268 C     31-37  cpm_*.dat (*=06,08...18)   Files are used as is
 269 C ----------------------------------------------------------------------------
 270 C
 271 C     DESCRIPTION OF INPUT PARAMETERS AND SWITCHES (FILE: hbt_parameters.in):
 272 C     ======================================================================
 273 C
 274 C     Control Switches:
 275 C     ================
 276 C
 277 C        ref_control       = 1 to read reference histograms from input files
 278 C                          = 2 to compute reference histograms by track
 279 C                              mixing from event pairs in the event input file.
 280 C
 281 C        switch_1d         = 0 to not compute the 1D two-body correlations.
 282 C                          = 1 to compute this using Q-invariant
 283 C                          = 2 to compute this using Q-total
 284 C                          = 3 to compute this using Q-3-vector
 285 C
 286 C        switch_3d         = 0 to not compute the 3D two-body correlations.
 287 C                          = 1 to compute this using the side-out-long form
 288 C                          = 2 to compute this using the YKP form.
 289 C
 290 C        switch_type       = 1 to fit only the like pair correlations
 291 C                          = 2 to fit only the unlike pair correlations
 292 C                          = 3 to fit both the like and unlike pair correl.
 293 C
 294 C        switch_coherence  = 0 for purely incoherent source (but can have
 295 C                              lambda < 1.0)
 296 C                          = 1 for mixed incoherent and coherent source
 297 C
 298 C        switch_coulomb    = 0 no Coulomb correction
 299 C                          = 1 Point source Gamow correction only
 300 C                          = 2 NA35 finite source size correction
 301 C                          = 3 Pratt/Cramer finite source size correction;
 302 C                              interpolated from input tables.
 303 C
 304 C        switch_fermi_bose =  1 Boson pairs
 305 C                          = -1 Fermion pairs
 306 C
 307 C        trk_accep         = 1.0 all adjusted tracks are written out
 308 C                          < 1.0 only this fraction, on average, of the
 309 C                            adjusted tracks are written out.  Used for
 310 C                            low multiplicity events.
 311 C
 312 C        print_full        = 0 for standard, minimal output
 313 C                          = 1 for full, comprehensive (large) output for
 314 C                              each event.
 315 C
 316 C        print_sector_data = 0 std. sector occupancy data printed out
 317 C                          = 1 to print sector occupancy and overflow info.
 318 C
 319 C     Particle ID and Search Parameters:
 320 C     =================================
 321 C
 322 C        n_pid_types       = 1 or 2 only, # particle types to correlate
 323 C
 324 C        pid(1), pid(2)    = Geant-3 particle ID code numbers
 325 C
 326 C        deltap            = maximum range for random momentum shifts in
 327 C                            GeV/c; px,py,pz independent; Def = 0.1 GeV/c.
 328 C
 329 C        maxit             = maximum # allowed iterations thru all the
 330 C                            tracks for each event. Default = 50.
 331 C                            If maxit=0, then calculate the correlations
 332 C                            for the input set of events without doing the
 333 C                            track adjustment procedure.
 334 C
 335 C        delchi            = min % change in total chi-square for which
 336 C                            the track adjustment iterations may stop,
 337 C                            Default = 0.1%.
 338 C
 339 C        irand             = initial random # seed, default = 12345
 340 C
 341 C     Source Function Parameters:
 342 C     ==========================
 343 C
 344 C        lambda            = Chaoticity
 345 C
 346 C        R_1d              = Spherical source model radius (fm)
 347 C
 348 C        Rside,Rout,Rlong  = Non-spherical Bertsch-Pratt source model (fm)
 349 C
 350 C        Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source
 351 C                            model (fm).
 352 C
 353 C        Q0                = NA35 Coulomb parameter for finite source size
 354 C                            in (GeV/c) - iff switch_coulomb = 2
 355 C                          = Spherical Coulomb source radius in (fm) iff
 356 C                            switch_coulomb = 3, used to interpolate the
 357 C                            input Pratt/Cramer discrete Coulomb source
 358 C                            radii tables.
 359 C
 360 C     One-body pT, phi, eta Acceptance Bins:
 361 C     =====================================
 362 C
 363 C        n_pt_bins, pt_min, pt_max  = # pt bins,  min/max pt accept. (GeV/c)
 364 C
 365 C        n_phi_bins,phi_min,phi_max = # phi bins, min/max phi accept. (deg.)
 366 C
 367 C        n_eta_bins,eta_min,eta_max = # eta bins, min/max eta accept.
 368 C
 369 C                                     [NOTE: For each the maximum # of bins
 370 C                                            must be .le. 100]
 371 C
 372 C     Two-body 1D and 3D Correlation Bins:
 373 C     ===================================
 374 C
 375 C        n_1d_fine,  binsize_1d_fine   = # and size (GeV/c), 1D - fine mesh
 376 C
 377 C        n_1d_coarse,binsize_1d_coarse = # and size (GeV/c), 1D - coarse mesh
 378 C
 379 C        n_3d_fine,  binsize_3d_fine   = # and size (GeV/c), 3D - fine mesh
 380 C
 381 C        n_3d_coarse,binsize_3d_coarse = # and size (GeV/c), 3D - coarse mesh
 382 C
 383 C                                      [NOTE: The maximum # of 1D bins (fine
 384 C                                             + coarse) must be .le. 100;
 385 C                                             The maximum # of 3D bins (either
 386 C                                             fine or coarse) must be .le.10).
 387 C                                             For both 1D and 3D there must be
 388 C                                             at least 1 fine bin and 1 coarse
 389 C                                             bin.]
 390 C        n_3d_fine_project             = # of 3D-fine bins to integrate over
 391 C                                        to form 1D projections.  This value
 392 C                                        must be .le. n_3d_fine.
 393 C
 394 C     Momentum Space Track-Sector Cells:
 395 C     =================================
 396 C
 397 C        n_px_bins,px_min,px_max = #, min,max px bins (GeV/c)
 398 C
 399 C        n_py_bins,py_min,py_max = #, min,max py bins (GeV/c)
 400 C
 401 C        n_pz_bins,pz_min,pz_max = #, min,max pz bins (GeV/c)
 402 C
 403 C                                  [NOTE: The maximum number of total sectors,
 404 C                                         equal to the product of the x-y-z
 405 C                                         number of cells must be .le.
 406 C                                         sec_maxlen which is defined in the
 407 C                                         /sec_trk_map/ data structure.]
 408 C
 409 C     Relative Chi-Square Weights:
 410 C     ===========================
 411 C
 412 C        chisq_wt_like_1d          = 1D, like pairs
 413 C        chisq_wt_unlike_1d        = 1D, unlike pairs
 414 C        chisq_wt_like_3d_fine     = 3D, like pairs, fine mesh
 415 C        chisq_wt_unlike_3d_fine   = 3D, unlike pairs, fine mesh
 416 C        chisq_wt_like_3d_coarse   = 3D, like pairs, coarse mesh
 417 C        chisq_wt_unlike_3d_coarse = 3D, unlike pairs, coarse mesh
 418 C        chisq_wt_hist1_1 = summed pt, phi, eta 1-body dist., PID#1
 419 C        chisq_wt_hist1_2 = summed pt, phi, eta 1-body dist., PID#2
 420 C
 421 C
 422 C     FORMAT for hbt_singles_reference.hist:
 423 C     =====================================
 424 C
 425 C     The output content for the one-body reference histograms is:
 426 C
 427 C     Line 1:  n_pid_types,pid(1),pid(2)
 428 C          2:  n_pt_bins,pt_min,pt_max
 429 C          3:  n_phi_bins,phi_min,phi_max
 430 C          4:  n_eta_bins,eta_min,eta_max
 431 C          5:  n_part_used_1_ref,n_part_used_2_ref
 432 C
 433 C     Then for PID #1:      (href1_pt_1(i),i=1,n_pt_bins)
 434 C     (One entry per line)  (href1_phi_1(i),i=1,n_phi_bins)
 435 C                           (href1_eta_1(i),i=1,n_eta_bins)
 436 C
 437 C     and for PID #2:       (href1_pt_2(i),i=1,n_pt_bins)
 438 C     (One entry per line)  (href1_phi_2(i),i=1,n_phi_bins)
 439 C                           (href1_eta_2(i),i=1,n_eta_bins)
 440 C
 441 C
 442 C     FORMAT for hbt_pair_reference.hist:
 443 C     ==================================
 444 C
 445 C     The output content for the two-body reference histograms is:
 446 C
 447 C     Line 1:  n_pid_types,pid(1),pid(2)
 448 C          2:  n_pt_bins,pt_min,pt_max
 449 C          3:  n_phi_bins,phi_min,phi_max
 450 C          4:  n_eta_bins,eta_min,eta_max
 451 C          5:  switch_1d,switch_3d,switch_type
 452 C          6:  n_1d_fine,n_1d_coarse,n_3d_fine,n_3d_coarse
 453 C          7:  binsize_1d_fine,binsize_1d_coarse,
 454 C              binsize_3d_fine,binsize_3d_coarse
 455 C          8:  num_pairs_like_ref,num_pairs_unlike_ref
 456 C
 457 C   The pair distributions (with one entry per line) are:
 458 C
 459 C   1D, like pairs:             (href_like_1d(i),i=1,n_1d_total)
 460 C
 461 C   1D, unlike pairs:           (href_unlike_1d(i),i=1,n_1d_total)
 462 C
 463 C   3D, like pairs, fine mesh:   href_like_3d_fine(i,j,k) ; (i,(j,(k,...)))
 464 C
 465 C   3D, like pairs, coarse mesh: href_like_3d_coarse(i,j,k) ; (i,(j,(k,...)))
 466 C
 467 C   3D, unlike, fine mesh:       href_unlike_3d_fine(i,j,k) ; (i,(j,(k,...)))
 468 C
 469 C   3D, unlike, coarse mesh:     href_unlike_3d_coarse(i,j,k) ; (i,(j,(k,...)))
 470 C
 471 C*************************************************************************
 472 C*************************************************************************
 473
 474       SUBROUTINE CTEST
 475       implicit none
 476
 477       Include 'common_parameters.inc'
 478       Include 'common_mesh.inc'
 479       Include 'common_histograms.inc'
 480       Include 'common_correlations.inc'
 481       Include 'common_coulomb.inc'
 482
 483       Include 'common_track.inc'
 484       Include 'common_track2.inc'
 485       Include 'common_sec_track.inc'
 486       Include 'common_sec_track2.inc'
 487       Include 'common_particle.inc'
 488
 489       write(*,*) ' '
 490       write(*,*) ' '
 491       write(*,*) ' '
 492
 493       write(*,*) 'Input data in Fort'
 494       write(*,*) ' '
 495       write(*,*) ' '
 496       write(*,*) ' '
 497       write(*,*) '      PARAMETERS'
 498       write(*,*) ' '
 499       write(*,*) 'ref_control',ref_control
 500       write(*,*) 'switch_1d',switch_1d
 501       write(*,*) 'switch_3d',switch_3d
 502       write(*,*) 'switch_type',switch_type
 503       write(*,*) 'switch_coherence',switch_coherence
 504       write(*,*) 'switch_coulomb',switch_coulomb
 505       write(*,*) 'switch_fermi_bose',switch_fermi_bose
 506       write(*,*) 'trk_accep',trk_accep
 507       write(*,*) 'print_full',print_full
 508       write(*,*) 'print_sector_data',print_sector_data
 509       write(*,*) 'n_pid_types',n_pid_types
 510       write(*,*) 'pid(1)', pid(1)
 511       write(*,*) 'pid(2)', pid(2)
 512       write(*,*) 'maxit',maxit
 513       write(*,*) 'irand',irand
 514       write(*,*) 'n_part_1_trk', n_part_1_trk
 515       write(*,*) 'n_part_2_trk ', n_part_2_trk
 516       write(*,*) 'n_part_tot_trk ', n_part_tot_trk
 517       write(*,*) 'n_part_used_1_trk ', n_part_used_1_trk
 518       write(*,*) 'n_part_used_2_trk',  n_part_used_2_trk
 519       write(*,*) 'n_part_1_trk2',    n_part_1_trk2
 520       write(*,*) 'n_part_2_trk2', n_part_2_trk2
 521       write(*,*) 'n_part_tot_trk2',   n_part_tot_trk2
 522       write(*,*) 'n_part_used_1_trk2',  n_part_used_1_trk2
 523       write(*,*) 'n_part_used_2_trk2',  n_part_used_2_trk2
 524       write(*,*) 'n_part_used_1_ref', n_part_used_1_ref
 525       write(*,*) 'n_part_used_2_ref ',  n_part_used_2_ref
 526       write(*,*) 'n_part_used_1_inc',  n_part_used_1_inc
 527       write(*,*) 'n_part_used_2_inc',     n_part_used_2_inc
 528       write(*,*) 'num_pairs_like',   num_pairs_like
 529       write(*,*) 'num_pairs_unlike',  num_pairs_unlike
 530       write(*,*) 'num_pairs_like_ref ', num_pairs_like_ref
 531       write(*,*) 'num_pairs_like_inc ',  num_pairs_like_inc
 532       write(*,*) 'num_pairs_unlike_inc', num_pairs_unlike_inc
 533       write(*,*) 'event_line_counter', event_line_counter
 534       write(*,*) 'file10_line_counter ',file10_line_counter
 535       write(*,*) 'lambda',lambda
 536       write(*,*) 'R_1d ',R_1d
 537       write(*,*) 'Rside',Rside
 538       write(*,*) 'Rout  ',  Rout
 539       write(*,*) 'Rlong  ', Rlong
 540       write(*,*) 'Rperp  ',  Rperp
 541       write(*,*) 'Rparallel  ', Rparallel
 542       write(*,*) 'R0 ',  R0
 543       write(*,*) 'Q0 ',  Q0
 544       write(*,*) 'deltap',deltap
 545       write(*,*) 'delchi',delchi
 546       write(*,*) 'pi ', pi
 547       write(*,*) 'rad ', rad
 548       write(*,*) 'hbc ', hbc
 549       write(*,*) 'chisq_wt_like_1d ', chisq_wt_like_1d
 550       write(*,*) 'chisq_wt_unlike_1d ',  chisq_wt_unlike_1d
 551       write(*,*) 'chisq_wt_like_3d_fine  ',chisq_wt_like_3d_fine
 552       write(*,*) 'chisq_wt_unlike_3d_fine ',  chisq_wt_unlike_3d_fine
 553       write(*,*) 'chisq_wt_like_3d_coarse ', chisq_wt_like_3d_coarse
 554       write(*,*) 'chisq_wt_unlike_3d_coarse',chisq_wt_unlike_3d_coarse
 555       write(*,*) 'chisq_wt_hist1_1 ', chisq_wt_hist1_1
 556       write(*,*) 'chisq_wt_hist1_2 ', chisq_wt_hist1_2
 557       write(*,*) ' '
 558       write(*,*) ' '
 559       write(*,*) ' '
 560       write(*,*) '         MESH       '
 561       write(*,*) ' '
 562       write(*,*) ' n_pt_bins ',  n_pt_bins
 563       write(*,*) ' pt_min ', pt_min
 564       write(*,*) ' pt_max  ', pt_max
 565       write(*,*) ' n_phi_bins ', n_phi_bins
 566       write(*,*) ' phi_min ', phi_min
 567       write(*,*) ' phi_max ', phi_max
 568       write(*,*) ' n_eta_bins ', n_eta_bins
 569       write(*,*) ' eta_min', eta_min
 570       write(*,*) ' eta_max', eta_max
 571       write(*,*) ' n_1d_fine ', n_1d_fine
 572       write(*,*) ' binsize_1d_fine ',binsize_1d_fine
 573       write(*,*) ' n_1d_coarse ',n_1d_coarse
 574       write(*,*) ' binsize_1d_coarse ', binsize_1d_coarse
 575       write(*,*) ' n_3d_fine  ',n_3d_fine
 576       write(*,*) ' binsize_3d_fine ',binsize_3d_fine
 577       write(*,*) ' n_3d_coarse ', n_3d_coarse
 578       write(*,*) ' binsize_3d_coarse ', binsize_3d_coarse
 579       write(*,*) ' n_3d_fine_project ',n_3d_fine_project
 580       write(*,*) ' n_px_bins ',n_px_bins
 581       write(*,*) ' px_min ',px_min
 582       write(*,*) ' px_max ', px_max
 583       write(*,*) ' n_py_bins ', n_py_bins
 584       write(*,*) ' py_min ',  py_min
 585       write(*,*) ' py_max',  py_max
 586       write(*,*) ' n_pz_bins ',n_pz_bins
 587       write(*,*) ' pz_min ', pz_min
 588       write(*,*) ' pz_max ', pz_max
 589       write(*,*) '  '
 590
 591       End
 592
 593
 594
 595       SUBROUTINE HBTPROCESSOR
 596       implicit none
 597
 598
 599       Include 'common_parameters.inc'
 600       Include 'common_mesh.inc'
 601       Include 'common_histograms.inc'
 602       Include 'common_correlations.inc'
 603       Include 'common_coulomb.inc'
 604
 605       Include 'common_track.inc'
 606       Include 'common_track2.inc'
 607       Include 'common_sec_track.inc'
 608       Include 'common_sec_track2.inc'
 609       Include 'common_particle.inc'
 610
 611 CCC   Set Data I/O control for ALICE or Standalone application
 612 C     ALICE = 1   ! This is for the ALICE AliRoot application
 613 CCC   ALICE = 0   ! This is for the standalone application
 614
 615 CCC   Initialize error code for ALICE application:
 616       errorcode = 0
 617
 618 CCC   Open Output Files:
 619
 620       open(unit=7,status='unknown',access='sequential',
 621      1     name='hbt_log.out')
 622       open(unit=8,status='unknown',access='sequential',
 623      1     name='hbt_simulation.out')
 624
 625 CCC   Initialize Arrays and Data Structures:
 626       If(ALICE .eq. 1) then
 627 C     In fact we not need to call initialization,
 628 C     because we can easily assume that is already done
 629          Call alihbtp_initialize
 630       Else If (ALICE .eq. 0) Then
 631          Call initialize
 632       End If
 633
 634       Write(6,100)
 635 CCC   Read Input Controls and Parameters:
 636       Call read_data(1)
 637       If(errorcode .eq. 1) Return
 638
 639 CCC   Setup values and check input parameter ranges and consistency:
 640       Call set_values
 641       If(errorcode .eq. 1) Return
 642
 643 CCC   Produce Basic Output File Header:
 644       Call write_data(1,0)
 645       If(errorcode .eq. 1) Return
 646
 647
 648       Write(6,101)
 649 CCC   Read Event Input file and fill flag files:
 650       Call read_data(2)
 651       If(errorcode .eq. 1) Return
 652
 653       Write(6,102)
 654 CCC   Get the Reference Histograms and write out if new calculation:
 655       Call getref_hist
 656       If(errorcode .eq. 1) Return
 657       Call write_data(3,0)
 658       If(errorcode .eq. 1) Return
 659       Write(6,103)
 660
 661       Write(6,104)
 662 CCC   Compute the correlation model and print out:
 663       Call correl_model
 664       Call write_data(4,0)
 665       If(errorcode .eq. 1) Return
 666
 667       Write(6,105)
 668 CCC   Carry out the Track Adjustment/Correlation Fitting Procedure:
 669       Call correlation_fit
 670       Write(6,106)
 671
 672 CCC   Final Output of Inclusive Quantities:
 673       Call write_data(6,0)
 674       If(errorcode .eq. 1) Return
 675
 676 CCC   Close Output Files:
 677       close(unit=7)
 678       close(unit=8)
 679
 680 CCC   Formats:
 681 100   Format(5x,'Read Input Controls, Setup values, check input:')
 682 101   Format(5x,'Read Event Input file and fill flag files:')
 683 102   Format(5x,'Get the Reference Histograms:')
 684 103   Format(5x,'Finished with Reference Histograms:')
 685 104   Format(5x,'Compute the correlation model:')
 686 105   Format(5x,'Start Track Adjustment/Correlation Fitting Procedure:')
 687 106   Format(5x,'Finished with Track Fitting Procedure:')
 688
 689       Return
 690       END
 691
 692 C-------------------------------------------------------------------
 693
 694
 695       subroutine initialize
 696       implicit none
 697
 698 CCC   This subroutine sets all arrays and structures to zero:
 699
 700       Include 'common_mesh.inc'
 701       Include 'common_histograms.inc'
 702       Include 'common_correlations.inc'
 703       Include 'common_coulomb.inc'
 704       Include 'common_event_summary.inc'
 705
 706       Include 'common_track.inc'
 707       Include 'common_track2.inc'
 708       Include 'common_sec_track.inc'
 709       Include 'common_sec_track2.inc'
 710       Include 'common_particle.inc'
 711
 712 CCC   Local Variable Type Declarations:
 713
 714       integer*4 i,j,k
 715
 716       do i = 1,trk_maxlen
 717          trk_id(i)            = 0
 718          trk_px_sec(i)        = 0
 719          trk_py_sec(i)        = 0
 720          trk_pz_sec(i)        = 0
 721          trk_sector(i)        = 0
 722          trk_flag(i)          = 0
 723          trk_out_flag(i)      = 0
 724          trk_merge_flag(i)    = 0
 725          trk_ge_pid(i)        = 0
 726          trk_start_vertex(i)  = 0
 727          trk_stop_vertex(i)   = 0
 728          trk_event_line(i)    = 0
 729          trk_px(i)            = 0.0
 730          trk_py(i)            = 0.0
 731          trk_pz(i)            = 0.0
 732          trk_E(i)             = 0.0
 733          trk_pt(i)            = 0.0
 734          trk_phi(i)           = 0.0
 735          trk_eta(i)           = 0.0
 736       end do
 737
 738       do i = 1,trk2_maxlen
 739          trk2_id(i)            = 0
 740          trk2_px_sec(i)        = 0
 741          trk2_py_sec(i)        = 0
 742          trk2_pz_sec(i)        = 0
 743          trk2_sector(i)        = 0
 744          trk2_flag(i)          = 0
 745          trk2_out_flag(i)      = 0
 746          trk2_merge_flag(i)    = 0
 747          trk2_ge_pid(i)        = 0
 748          trk2_start_vertex(i)  = 0
 749          trk2_stop_vertex(i)   = 0
 750          trk2_event_line(i)    = 0
 751          trk2_px(i)            = 0.0
 752          trk2_py(i)            = 0.0
 753          trk2_pz(i)            = 0.0
 754          trk2_E(i)             = 0.0
 755          trk2_pt(i)            = 0.0
 756          trk2_phi(i)           = 0.0
 757          trk2_eta(i)           = 0.0
 758       end do
 759
 760       do i = 1,sec_maxlen
 761          stm_sec_id(i)         = 0
 762          stm_n_trk_sec(i)      = 0
 763          stm_flag(i)           = 0
 764          do j = 1,max_trk_sec
 765             stm_track_id(j,i)  = 0
 766          end do
 767       end do
 768
 769       do i = 1,sec_maxlen2
 770          stm2_sec_id(i)         = 0
 771          stm2_n_trk_sec(i)      = 0
 772          stm2_flag(i)           = 0
 773          do j = 1,max_trk_sec2
 774             stm2_track_id(j,i)  = 0
 775          end do
 776       end do
 777
 778       do i = 1,part_maxlen
 779          part_id(i)       = 0
 780          part_charge(i)   = 0
 781          part_mass(i)     = 0.0
 782          part_lifetime(i) = 0.0
 783       end do
 784
 785       do i = 1,max_trk_save
 786          old_sec_trkid(i) = 0
 787          new_sec_trkid(i) = 0
 788       end do
 789
 790       do i = 1,max_h_1d
 791          hist_like_1d(i)   = 0
 792          hist_unlike_1d(i) = 0
 793          htmp_like_1d(i)   = 0
 794          htmp_unlike_1d(i) = 0
 795          href_like_1d(i)   = 0
 796          href_unlike_1d(i) = 0
 797          hinc_like_1d(i)   = 0
 798          hinc_unlike_1d(i) = 0
 799          hist1_pt_1(i)     = 0
 800          hist1_pt_2(i)     = 0
 801          hist1_phi_1(i)    = 0
 802          hist1_phi_2(i)    = 0
 803          hist1_eta_1(i)    = 0
 804          hist1_eta_2(i)    = 0
 805          htmp1_pt_1(i)     = 0
 806          htmp1_pt_2(i)     = 0
 807          htmp1_phi_1(i)    = 0
 808          htmp1_phi_2(i)    = 0
 809          htmp1_eta_1(i)    = 0
 810          htmp1_eta_2(i)    = 0
 811          href1_pt_1(i)     = 0
 812          href1_pt_2(i)     = 0
 813          href1_phi_1(i)    = 0
 814          href1_phi_2(i)    = 0
 815          href1_eta_1(i)    = 0
 816          href1_eta_2(i)    = 0
 817          hinc1_pt_1(i)     = 0
 818          hinc1_pt_2(i)     = 0
 819          hinc1_phi_1(i)    = 0
 820          hinc1_phi_2(i)    = 0
 821          hinc1_eta_1(i)    = 0
 822          hinc1_eta_2(i)    = 0
 823       end do
 824
 825       do i = 1,max_h_3d
 826       do j = 1,max_h_3d
 827       do k = 1,max_h_3d
 828          hist_like_3d_fine(i,j,k)     = 0
 829          hist_unlike_3d_fine(i,j,k)   = 0
 830          hist_like_3d_coarse(i,j,k)   = 0
 831          hist_unlike_3d_coarse(i,j,k) = 0
 832          htmp_like_3d_fine(i,j,k)     = 0
 833          htmp_unlike_3d_fine(i,j,k)   = 0
 834          htmp_like_3d_coarse(i,j,k)   = 0
 835          htmp_unlike_3d_coarse(i,j,k) = 0
 836          href_like_3d_fine(i,j,k)     = 0
 837          href_unlike_3d_fine(i,j,k)   = 0
 838          href_like_3d_coarse(i,j,k)   = 0
 839          href_unlike_3d_coarse(i,j,k) = 0
 840          hinc_like_3d_fine(i,j,k)     = 0
 841          hinc_unlike_3d_fine(i,j,k)   = 0
 842          hinc_like_3d_coarse(i,j,k)   = 0
 843          hinc_unlike_3d_coarse(i,j,k) = 0
 844       end do
 845       end do
 846       end do
 847
 848       do i = 1,max_c2_1d
 849          c2mod_like_1d(i)   = 0.0
 850          c2mod_unlike_1d(i) = 0.0
 851          c2fit_like_1d(i)   = 0.0
 852          c2fit_unlike_1d(i) = 0.0
 853          c2err_like_1d(i)   = 0.0
 854          c2err_unlike_1d(i) = 0.0
 855       end do
 856
 857       do i = 1,max_c2_3d
 858       do j = 1,max_c2_3d
 859       do k = 1,max_c2_3d
 860          c2mod_like_3d_fine(i,j,k)     = 0.0
 861          c2mod_unlike_3d_fine(i,j,k)   = 0.0
 862          c2mod_like_3d_coarse(i,j,k)   = 0.0
 863          c2mod_unlike_3d_coarse(i,j,k) = 0.0
 864          c2fit_like_3d_fine(i,j,k)     = 0.0
 865          c2fit_unlike_3d_fine(i,j,k)   = 0.0
 866          c2fit_like_3d_coarse(i,j,k)   = 0.0
 867          c2fit_unlike_3d_coarse(i,j,k) = 0.0
 868          c2err_like_3d_fine(i,j,k)     = 0.0
 869          c2err_unlike_3d_fine(i,j,k)   = 0.0
 870          c2err_like_3d_coarse(i,j,k)   = 0.0
 871          c2err_unlike_3d_coarse(i,j,k) = 0.0
 872       end do
 873       end do
 874       end do
 875
 876       do i = 1,max_c2_coul
 877          c2_coul_like(i)   = 0.0
 878          c2_coul_unlike(i) = 0.0
 879          q_coul(i)         = 0.0
 880       end do
 881
 882       do i = 1,max_events
 883          num_iter(i)                      = 0.0
 884          n_part_used_1_store(i)           = 0.0
 885          n_part_used_2_store(i)           = 0.0
 886          num_sec_flagged_store(i)         = 0.0
 887          frac_trks_out(i)                 = 0.0
 888          frac_trks_flag(i)                = 0.0
 889          chisq_like_1d_store(i)           = 0.0
 890          chisq_unlike_1d_store(i)         = 0.0
 891          chisq_like_3d_fine_store(i)      = 0.0
 892          chisq_unlike_3d_fine_store(i)    = 0.0
 893          chisq_like_3d_coarse_store(i)    = 0.0
 894          chisq_unlike_3d_coarse_store(i)  = 0.0
 895          chisq_hist1_1_store(i)           = 0.0
 896          chisq_hist1_2_store(i)           = 0.0
 897          chisq_total_store(i)             = 0.0
 898       end do
 899
 900       Return
 901       END
 902
 903 C---------------------------------------------------------------------
 904
 905
 906       subroutine set_values
 907       implicit none
 908
 909 CCC   This subroutine sets parameters based on the main input.
 910 CCC   The consistency of the input parameters and controls is
 911 CCC   checked.  Any problems are reported in the Log File,
 912 CCC   'hbt_log.out'.  Most inconsistencies or array size limit
 913 CCC   overflows will cause the code execution to STOP.
 914
 915       Include 'common_parameters.inc'
 916       Include 'common_mesh.inc'
 917       Include 'common_histograms.inc'
 918       Include 'common_correlations.inc'
 919       Include 'common_coulomb.inc'
 920
 921       Include 'common_track.inc'
 922       Include 'common_track2.inc'
 923       Include 'common_sec_track.inc'
 924       Include 'common_sec_track2.inc'
 925       Include 'common_particle.inc'
 926
 927 CCC   Local Variable Type Declarations:
 928
 929       integer*4 iphistep, ptmaxsteps, iptstep
 930
 931       real*4 px1,py1,pz1,E1, pt1,phi1
 932       real*4 px2,py2,pz2,E2
 933       real*4 pt_step,phi_step
 934       real*4 pxstepmin, pxstepmax, pystepmin, pystepmax
 935
 936 CCC  Check Input Controls:
 937
 938       if(ref_control .lt. 1 .or. ref_control .gt. 2) then
 939          write(7,101) ref_control
 940          errorcode = 1
 941          Return
 942       end if
 943
 944       if(switch_1d .lt. 0 .or. switch_1d .gt. 3) then
 945          write(7,102) switch_1d
 946          errorcode = 1
 947          Return
 948       end if
 949
 950       if(switch_3d .lt. 0 .or. switch_3d .gt. 2) then
 951          write(7,103) switch_3d
 952          errorcode = 1
 953          Return
 954       end if
 955
 956       if(switch_type .lt. 1 .or. switch_type .gt. 3) then
 957          write(7,104) switch_type
 958          errorcode = 1
 959          Return
 960       end if
 961
 962       if(switch_coherence .lt. 0 .or. switch_coherence .gt. 1) then
 963          write(7,105) switch_coherence
 964          errorcode = 1
 965          Return
 966       end if
 967
 968       if(switch_coulomb .lt. 0 .or. switch_coulomb .gt. 3) then
 969          write(7,106) switch_coulomb
 970          errorcode = 1
 971          Return
 972       end if
 973
 974       if(switch_fermi_bose.ne.-1 .and. switch_fermi_bose.ne.1) then
 975          write(7,107) switch_fermi_bose
 976          errorcode = 1
 977          Return
 978       end if
 979
 980       if(n_pid_types .lt. 1 .or. n_pid_types .gt. 2) then
 981          write(7,108) n_pid_types
 982          errorcode = 1
 983          Return
 984       end if
 985
 986       if(switch_type .ge. 2 .and. n_pid_types .eq. 1) then
 987          write(7,109) switch_type, n_pid_types
 988          errorcode = 1
 989          Return
 990       end if
 991
 992       if(n_pid_types .eq. 1) then
 993          if(pid(1).gt.0 .and. pid(2).gt.0) then
 994             write(7,1091) pid(1),pid(2)
 995             errorcode = 1
 996             Return
 997          end if
 998       end if
 999
1000       if(pid(1).eq.0 .and. pid(2).eq.0) then
1001          write(7,1092)
1002          errorcode = 1
1003          Return
1004       end if
1005
1006       if(n_pid_types .eq. 2) then
1007          if(pid(1).gt.0.and.pid(2).gt.0.and.pid(1).ne.pid(2))then
1008             continue
1009          else
1010             write(7,1093) pid(1), pid(2)
1011             errorcode = 1
1012             Return
1013          end if
1014       end if
1015
1016       if(pid(1).gt.0.and.pid(2).gt.0.and.pid(1).eq.pid(2))then
1017          write(7,1094) pid(1), pid(2)
1018          errorcode = 1
1019          Return
1020       end if
1021
1022       if(trk_accep .le. 0.0) then
1023          write(7,10941) trk_accep
1024          errorcode = 1
1025          Return
1026
1027       end if
1028
1029 CCC   Check Input Parameters:
1030
1031       if(deltap .le. 0.0) deltap = 0.1
1032       if(maxit  .lt. 0  ) maxit = 50
1033       if(delchi .lt. 0.0) delchi = 0.1
1034       if(irand  .le. 0  ) irand = 12345
1035
1036 CCC   Check Coulomb source radius in range for Pratt type Coulomb correction.
1037
1038       if(switch_coulomb .eq. 3 .and. (Q0 .lt. coulradmin .or.
1039      1   Q0 .gt. coulradmax)) then
1040          write(7,132) Q0
1041          errorcode = 1
1042          Return
1043       end if
1044
1045 CCC   Load the Pratt type Coulomb correction if this form is selected.
1046
1047       if(switch_coulomb .eq. 3 .and. (Q0 .ge. coulradmin .and.
1048      1   Q0 .le. coulradmax)) then
1049          Call read_data(6)
1050       end if
1051
1052 CCC   Check and determine the one-body distribution's binning:
1053
1054       if(n_pt_bins .lt. 1 .or. n_pt_bins .gt. max_h_1d) then
1055          write(7,110) n_pt_bins
1056          errorcode = 1
1057          Return
1058       end if
1059
1060       if(n_phi_bins .lt. 1 .or. n_phi_bins .gt. max_h_1d) then
1061          write(7,111) n_phi_bins
1062          errorcode = 1
1063          Return
1064       end if
1065
1066       if(n_eta_bins .lt. 1 .or. n_eta_bins .gt. max_h_1d) then
1067          write(7,112) n_eta_bins
1068          errorcode = 1
1069          Return
1070       end if
1071
1072       if(pt_min .gt. pt_max .or. pt_min .lt. 0.0) then
1073          write(7,113) pt_min, pt_max
1074          errorcode = 1
1075          Return
1076       end if
1077
1078       if(phi_min.gt.phi_max .or. phi_min.lt.0.0 .or.
1079      1   phi_max.gt.360.0) then
1080          write(7,114) phi_min, phi_max
1081          errorcode = 1
1082          Return
1083       end if
1084
1085       if(eta_min .gt. eta_max) then
1086          write(7,115) eta_min, eta_max
1087          errorcode = 1
1088          Return
1089       end if
1090
1091       pt_bin_size  = (pt_max  - pt_min )/float(n_pt_bins)
1092       phi_bin_size = (phi_max - phi_min)/float(n_phi_bins)
1093       eta_bin_size = (eta_max - eta_min)/float(n_eta_bins)
1094
1095 CCC   Check and determine the two-body distribution's binning:
1096
1097       n_1d_total = n_1d_fine + n_1d_coarse
1098       n_3d_total = n_3d_fine + n_3d_coarse - 1
1099
1100       if(switch_1d .gt. 0) then
1101       if(n_1d_fine .lt. 1) then
1102          write(7,116) n_1d_fine
1103          errorcode = 1
1104          Return
1105       end if
1106
1107       if(n_1d_coarse .lt. 1) then
1108          write(7,117) n_1d_coarse
1109          errorcode = 1
1110          Return
1111       end if
1112
1113       if(n_1d_total .gt. max_h_1d) then
1114          write(7,118) n_1d_total
1115          errorcode = 1
1116          Return
1117       end if
1118
1119       qmid_1d = binsize_1d_fine  *float(n_1d_fine)
1120       qmax_1d = binsize_1d_coarse*float(n_1d_coarse) + qmid_1d
1121       end if
1122
1123       if(switch_3d .gt. 0) then
1124       if(n_3d_fine .lt. 1 .or. n_3d_fine .gt. max_h_3d) then
1125          write(7,119) n_3d_fine
1126          errorcode = 1
1127          Return
1128       end if
1129
1130       if(n_3d_coarse .lt. 1 .or. n_3d_coarse .gt. max_h_3d) then
1131          write(7,120) n_3d_coarse
1132          errorcode = 1
1133          Return
1134       end if
1135
1136       qmid_3d = binsize_3d_fine  *float(n_3d_fine)
1137       qmax_3d = binsize_3d_coarse*float(n_3d_coarse)
1138
1139       if(abs(qmid_3d - binsize_3d_coarse) .gt. 0.00001) then
1140          write(7,121) qmid_3d, binsize_3d_coarse
1141          errorcode = 1
1142          Return
1143       end if
1144
1145       if(n_3d_fine_project .gt. n_3d_fine) then
1146          write(7,1211) n_3d_fine_project, n_3d_fine
1147          n_3d_fine_project = n_3d_fine
1148       end if
1149
1150       if(n_3d_fine_project .lt. 1) then
1151          write(7,1212) n_3d_fine_project
1152          n_3d_fine_project = 1
1153       end if
1154       end if
1155
1156 CCC   Check and determine Track-Sector Parameters:
1157
1158       if(n_px_bins .lt. 1) then
1159          write(7,122) n_px_bins
1160          errorcode = 1
1161          Return
1162       end if
1163
1164       if(n_py_bins .lt. 1) then
1165          write(7,123) n_py_bins
1166          errorcode = 1
1167          Return
1168       end if
1169
1170       if(n_pz_bins .lt. 1) then
1171          write(7,124) n_pz_bins
1172          errorcode = 1
1173          Return
1174       end if
1175
1176       n_sectors = n_px_bins * n_py_bins * n_pz_bins
1177       if(n_sectors .gt. sec_maxlen) then
1178          write(7,125) n_sectors
1179          errorcode = 1
1180          Return
1181       end if
1182
1183       if(n_sectors .gt. sec_maxlen2 .and. ref_control .eq. 2) then
1184          write(7,1251) n_sectors
1185          errorcode = 1
1186          Return
1187       end if
1188
1189       if(trk_maxlen .ne. trk2_maxlen .and. ref_control .eq. 2) then
1190          write(7,1252)
1191          errorcode = 1
1192          Return
1193       end if
1194
1195       if(max_trk_save .ne. max_trk_sec .or.
1196      1   max_trk_save .ne. max_trk_sec2 .or.
1197      2   max_trk_sec  .ne. max_trk_sec2) then
1198          write(7,12521) max_trk_save,max_trk_sec,max_trk_sec2
1199          errorcode = 1
1200          Return
1201       end if
1202
1203       delpx = (px_max - px_min)/float(n_px_bins)
1204       delpy = (py_max - py_min)/float(n_py_bins)
1205       delpz = (pz_max - pz_min)/float(n_pz_bins)
1206
1207 CCC   Check that the Track-Sector Grid includes the acceptance range:
1208 CCC   The Track-Sector Grid is a 3D {px,py,pz} box, while the acceptance
1209 CCC   is defined in cylindrical {pt,phi,eta} coordinates.
1210 CCC
1211 CCC   Check the z-momentum components:
1212
1213       if(eta_min .ge. 0.0) then
1214          Call Hbtp_kin(px1,py1,pz1,E1,pt_min,0.0,eta_min,0.14,2)
1215          Call Hbtp_kin(px2,py2,pz2,E2,pt_max,0.0,eta_max,0.14,2)
1216       else if(eta_max .le. 0.0) then
1217          Call Hbtp_kin(px1,py1,pz1,E1,pt_max,0.0,eta_min,0.14,2)
1218          Call Hbtp_kin(px2,py2,pz2,E2,pt_min,0.0,eta_max,0.14,2)
1219       else if(eta_min .le. 0.0 .and. eta_max .ge. 0.0) then
1220          Call Hbtp_kin(px1,py1,pz1,E1,pt_max,0.0,eta_min,0.14,2)
1221          Call Hbtp_kin(px2,py2,pz2,E2,pt_max,0.0,eta_max,0.14,2)
1222       end if
1223
1224       if(pz1 .lt. pz_min .or. pz2 .gt. pz_max) then
1225          write(7,126) pz1,pz_min,pz2,pz_max
1226          errorcode = 1
1227          Return
1228       end if
1229
1230 CCC   Check the x,y-momentum components by scanning over the perimeter
1231 CCC   of the (pt,phi) acceptance domain space with 100 trial grid points.
1232 CCC   The overall required px_min, px_max, py_min and py_max to cover the
1233 CCC   acceptance by the track-sectors is determined.  These values are
1234 CCC   then compared with the min/max px and py ranges for the track-sectors.
1235 CCC
1236 CCC   Divide the pt and phi acceptance ranges into 24 equal steps:
1237
1238       pt_step  = (pt_max  - pt_min)/24.0
1239       phi_step = (phi_max - phi_min)/24.0
1240       pxstepmax = -1000.
1241       pxstepmin =  1000.
1242       pystepmax = -1000.
1243       pystepmin =  1000.
1244       phi1 = phi_min - phi_step
1245       do iphistep = 1,25
1246          phi1 = phi1 + phi_step
1247          ptmaxsteps = 2
1248          if(iphistep.eq.1 .or. iphistep.eq.25) ptmaxsteps = 25
1249          pt1 = pt_min - pt_step
1250          do iptstep = 1,ptmaxsteps
1251             if(iphistep.eq.1 .or. iphistep.eq.25) then
1252                pt1 = pt1 + pt_step
1253             else if(iphistep.gt.1 .and. iphistep.lt.25) then
1254                if(iptstep.eq.1) pt1 = pt_min
1255                if(iptstep.eq.2) pt1 = pt_max
1256             end if
1257             Call Hbtp_kin(px1,py1,pz1,E1,pt1,phi1,0.0,0.14,2)
1258                if(px1.gt.pxstepmax) pxstepmax = px1
1259                if(px1.lt.pxstepmin) pxstepmin = px1
1260                if(py1.gt.pystepmax) pystepmax = py1
1261                if(py1.lt.pystepmin) pystepmin = py1
1262          end do
1263       end do
1264
1265       if(pxstepmin .lt. px_min .or. pxstepmax .gt. px_max) then
1266          write(7,127) pxstepmin,px_min,pxstepmax,px_max
1267          errorcode = 1
1268          Return
1269       end if
1270
1271       if(pystepmin .lt. py_min .or. pystepmax .gt. py_max) then
1272          write(7,128) pystepmin,py_min,pystepmax,py_max
1273          errorcode = 1
1274          Return
1275       end if
1276
1277 CCC   Load Geant Particle Data:
1278       Call Hbtp_particle_prop
1279
1280 CCC   Check Requested Particle ID Numbers:
1281
1282       if(n_pid_types.eq.1 .and. pid(1).le.0 .and. pid(2).le.0) then
1283          write(7,131) pid(1),pid(2)
1284          errorcode = 1
1285          Return
1286       end if
1287
1288 CCC   Initialize Masses to 0.0
1289
1290       mass1 = 0.0
1291       mass2 = 0.0
1292
1293       if(n_pid_types .eq. 1 .and. pid(1) .ne. 0) then
1294          if(pid(1) .lt. 1 .or. pid(1) .gt. part_maxlen) then
1295             write(7,129) pid(1)
1296             errorcode = 1
1297             Return
1298          else
1299             mass1 = part_mass(pid(1))
1300          end if
1301       else if(n_pid_types .eq. 1 .and. pid(2) .ne. 0) then
1302          if(pid(2) .lt. 1 .or. pid(2) .gt. part_maxlen) then
1303             write(7,130) pid(2)
1304             errorcode = 1
1305             Return
1306          else
1307             mass2 = part_mass(pid(2))
1308          end if
1309       else if(n_pid_types .eq. 2) then
1310          if(pid(1) .lt. 1 .or. pid(1) .gt. part_maxlen) then
1311             write(7,129) pid(1)
1312             errorcode = 1
1313             Return
1314          else
1315             mass1 = part_mass(pid(1))
1316          end if
1317          if(pid(2) .lt. 1 .or. pid(2) .gt. part_maxlen) then
1318             write(7,130) pid(2)
1319             errorcode = 1
1320             Return
1321          else
1322             mass2 = part_mass(pid(2))
1323          end if
1324       end if
1325
1326 CCC   Set Math Constants:
1327
1328       pi = 3.141592654
1329       hbc = 0.19732891
1330       rad = 180.0/pi
1331
1332 CCC   FORMATS:
1333
1334 101   Format(5x,'ref_control = ',I5,'Out of Range - STOP')
1335 102   Format(5x,'switch_1d   = ',I5,'Out of Range - STOP')
1336 103   Format(5x,'switch_3d   = ',I5,'Out of Range - STOP')
1337 104   Format(5x,'switch_type = ',I5,'Out of Range - STOP')
1338 105   Format(5x,'switch_coherence = ',I5,'Out of Range - STOP')
1339 106   Format(5x,'switch_coulomb = ',I5,'Out of Range - STOP')
1340 107   Format(5x,'switch_fermi_bose = ',I5,'Out of Range - STOP')
1341 108   Format(5x,'n_pid_types = ',I5,'Out of Range - STOP')
1342 109   Format(5x,'switch_type & n_pid_types = ',2I5,
1343      1       'Incompatible - STOP')
1344 1091  Format(5x,'For n_pid_types=1, cannot have both PID#1,2 = ',
1345      1       2I5,' .ne.0 - STOP')
1346 1092  Format(5x,'Both PID# 1 and 2 = 0, - STOP')
1347 1093  Format(5x,'For n_pid_types=2, PID#1,2 = ',2I5,
1348      1      ' are incorrect - STOP')
1349 1094  Format(5x,'Both PID# 1,2 = ',2I5,' are equal - STOP')
1350 10941 Format(5x,'Track acceptance output frac .le. 0.0 - STOP')
1351 132   Format(5x,'Coulomb source radius = ',E12.4,' - For Pratt ',
1352      1       'Correction, Out of Range - STOP')
1353 110   Format(5x,'# pt bins  = ',I5,'Out of Range - STOP')
1354 111   Format(5x,'# phi bins = ',I5,'Out of Range - STOP')
1355 112   Format(5x,'# eta bins = ',I5,'Out of Range - STOP')
1356 113   Format(5x,'pt min/max accept. range = ',2E12.4,' is wrong-STOP')
1357 114   Format(5x,'phi min/max accept. range = ',2E12.4,' is wrong-STOP')
1358 115   Format(5x,'eta min/max accept. range = ',2E12.4,' is wrong-STOP')
1359 116   Format(5x,'# 1d fine mesh for C2 = ',I5,' .lt.1 - STOP')
1360 117   Format(5x,'# 1d coarse mesh for C2 = ',I5,' .lt.1 - STOP')
1361 118   Format(5x,'Total # 1d mesh for C2 = ',I5,' .gt.max_h_1d - STOP')
1362 119   Format(5x,'# 3d fine mesh for C2 = ',I5,'Out of Range - STOP')
1363 120   Format(5x,'# 3d coarse mesh for C2 = ',I5,'Out of Range - STOP')
1364 121   Format(5x,'3D C2 fine range & coarse grid = ',2E12.4,
1365      1       'Not Equal - STOP')
1366 1211  Format(5x,'# 3D fine bins projected = ',I5,
1367      1       ' TOO BIG - reduce to n_3d_fine = ',I5)
1368 1212  Format(5x,'# 3D fine bins projected = ',I5,
1369      1       ' Set to 1')
1370 122   Format(5x,'#track-sector px bins = ',I5,' .lt.1 - STOP')
1371 123   Format(5x,'#track-sector py bins = ',I5,' .lt.1 - STOP')
1372 124   Format(5x,'#track-sector pz bins = ',I5,' .lt.1 - STOP')
1373 125   Format(5x,'Total # trk-sec = ',I5,' .gt.sec_maxlen - STOP')
1374 1251  Format(5x,'Total # trk-sec = ',I5,' .gt.sec_maxlen2 for ',
1375      1       'Reference calc. - STOP')
1376 1252  Format(5x,'trk_maxlen .ne. trk2_maxlen for Ref. Calc. - STOP')
1377 12521 Format(5x,'max_trk_save,max_trk_sec,max_trk_sec2 = ',
1378      1       3I5,' are not all equal - STOP')
1379 126   Format(5x,'pz accept. not covered by sectors-STOP:',4E12.4)
1380 127   Format(5x,'px accept. not covered by sectors-STOP:',4E12.4)
1381 128   Format(5x,'py accept. not covered by sectors-STOP:',4E12.4)
1382 131   Format(5x,'Particle ID values = ',2I5,' not valid - STOP')
1383 129   Format(5x,'Particle ID value #1 = ',I5,' not valid - STOP')
1384 130   Format(5x,'Particle ID value #2 = ',I5,' not valid - STOP')
1385
1386       Return
1387       END
1388
1389 C----------------------------------------------------------------------
1390
1391
1392       subroutine read_data(mode)
1393       implicit none
1394
1395 CCC   This subroutine does all the data input associated with all input
1396 C     files.  Some diagnostic output is printed here if errors occur
1397 C     during the file reading.  Two auxiliary output files, which tag
1398 C     the events input tracks are written out.
1399 C
1400 C     The type of input is controlled by the value of 'mode'
1401 C     where:
1402 C     (The following mostly applies to the standalone application
1403 C     that reads from files and writes temporary scratch files.
1404 C     This is the ALICE=0 mode.)
1405 C
1406 C        mode = 1, read the parameter and switches input file
1407 C
1408 C        mode = 2, scan the event text file and write out two
1409 C                  auxiliary output/tag files; select and mark
1410 C                  accepted tracks to use.
1411 C
1412 C        mode = 3, read the reference pair and one-body histograms
1413 C
1414 C        mode = 4, read the next event from the event text file,
1415 C                  'event_text.in,' and load the accepted tracks
1416 C                  into the 'trk' data structure.
1417 C
1418 C        mode = 5, same as mode=4, except the accepted tracks are
1419 C                  loaded into the 'trk2' data structure.
1420 C
1421 C        mode = 6, read the input Coulomb correction tables and
1422 C                  interpolate for the requested source radius, arrays
1423 C                  in common/coulomb/ are filled for like and unlike
1424 C                  charged pairs.
1425 C
1426 C        mode = 7, read the next event from the event text file,
1427 C                  'event_text.in,' and load the accepted tracks
1428 C                  into the 'trk' data structure. Then copy the event
1429 C                  data in 'event_text.in' to 'event_text_aux.in' and
1430 C                  from 'event_tracks.select' to 'event_tracks_aux.select'
1431 C
1432 C        mode = 8, read contents of 'event_text_aux.in' using flag values
1433 C                  in 'event_tracks_aux.select' and copy into
1434 C                  'event_hbt_text.out' (i.e. the main event output file)
1435 C                  where the momentum values for accepted tracks are
1436 C                  overwritten with the adjusted (correlated) parameters
1437 C                  in the 'trk' data structure.
1438 C
1439 C                  If trk_accep .lt. 1.0, then only write this fraction
1440 C                  of the final tracks, as determined by a random number
1441 C                  throw.
1442 C
1443 C     Summary of Files:
1444 C     ----------------
1445 C
1446 C File Unit #        Filename               Description
1447 C ---------------------------------------------------------------------------
1448 C     1       hbt_parameters.in          Input switches, parameters
1449 C     2       event_text.in              Event Gen output, GSTAR text format
1450 C     3       event_line.flags           Event file line flags
1451 C     4       event_tracks.select        Event file selected tracks
1452 C     7       hbt_log.out                Log and error messages
1453 C     8       hbt_simulation.out         Full Output
1454 C     9       hbt_pair_reference.hist    Reference pair histograms
1455 C    10       event_hbt_text.out         Updated/correlated event text file
1456 C    11       hbt_singles_reference.hist Reference one-body histograms
1457 C    12       event_text_aux.in          Tmp. copy of 'event_text.in'/event
1458 C    14       event_tracks_aux.select    Tmp. copy 'event_tracks.select'/event
1459 C   21-27     cpp_*.dat (*=06,08,...18)  Like pair Pratt Coul. Correct
1460 C   31-37     cpm_*.dat (*=06,08,...18)  Unlike pair Pratt Coul. Correct
1461 C ---------------------------------------------------------------------------
1462 C
1463
1464       Include 'common_parameters.inc'
1465       Include 'common_mesh.inc'
1466       Include 'common_histograms.inc'
1467       Include 'common_correlations.inc'
1468       Include 'common_coulomb.inc'
1469
1470       Include 'common_track.inc'
1471       Include 'common_track2.inc'
1472       Include 'common_particle.inc'
1473
1474       integer LNBLNK
1475
1476 CCC   Local Variable Type Declarations:
1477
1478       real*4 px,py,pz,E,pt,phi,eta,mass
1479       real*4 acheck(10), function(20)
1480       real*4 hbtpran
1481
1482       integer*4 i,j,k,mode,flag,flag4,flag0,ntracks
1483       integer*4 ge_pid,tid,start_v,stop_v,eg_pid
1484       integer*4 ref_check,pidok,accepok,check(13)
1485       integer*4 event_counter,track_counter
1486       integer*4 track_counter_1,track_counter_2
1487
1488       character*5 evg_label,event_line,vertex_line,track_line,dummy
1489       character*5 gener_line
1490       character*80 comment_event_label
1491       character*87 vertex_label
1492       character*93 gener_label
1493
1494       parameter (event_line  = 'EVENT')
1495       parameter (vertex_line = 'VERTE')
1496       parameter (track_line  = 'TRACK')
1497       parameter (gener_line  = 'GENER')
1498       parameter (flag4 = 4)
1499       parameter (flag0 = 0)
1500 C     ALICE USE ONLY
1501       CHARACTER*100 CHROOT
1502       CHARACTER*100 FILNAM
1503       INTEGER*4 LNROOT
1504       LOGICAL EXISTS
1505       CHROOT=' '
1506 C
1507
1508 CCC   Begin Data Input Options:
1509
1510 C------------------------
1511       IF (mode.eq.1) THEN              ! Read Input parameters from File#1
1512 C------------------------
1513
1514 CCC   For standalone version (ALICE = 0), read parameters from
1515 CCC   File Unit=1, 'hbt_parameters.in'
1516 CCC   For ALICE-ROOT version (ALICE=1) load parameters from Call to C++ funct
1517       If(ALICE .eq. 1) then
1518          Call AliHbtp_SetParameters
1519       Else If(ALICE .eq. 0) Then
1520
1521          open(unit=1,type='old',access='sequential',
1522      1        name='hbt_parameters.in')
1523
1524 CCC   Read Control Switches:  (See Main program listing for complete
1525 CCC                            description of input parameters)
1526
1527       read(1,*) ref_control
1528       read(1,*) switch_1d
1529       read(1,*) switch_3d
1530       read(1,*) switch_type
1531       read(1,*) switch_coherence
1532       read(1,*) switch_coulomb
1533       read(1,*) switch_fermi_bose
1534       read(1,*) trk_accep
1535       read(1,*) print_full
1536       read(1,*) print_sector_data
1537
1538 CCC   Read Parameters:
1539
1540       read(1,*) n_pid_types
1541       read(1,*) pid(1),pid(2)
1542       read(1,*) deltap
1543       read(1,*) maxit
1544       read(1,*) delchi
1545       read(1,*) irand
1546
1547 CCC   Read Source Parameters:
1548
1549       read(1,*) lambda
1550       read(1,*) R_1d
1551       read(1,*) Rside, Rout, Rlong
1552       read(1,*) Rperp, Rparallel, R0
1553       read(1,*) Q0
1554
1555 CCC   Read one-body {pt,phi,eta} bins:
1556
1557       read(1,*) n_pt_bins ,pt_min ,pt_max
1558       read(1,*) n_phi_bins,phi_min,phi_max
1559       read(1,*) n_eta_bins,eta_min,eta_max
1560
1561 CCC   Read two-body 1D and 3D bins:
1562
1563       read(1,*) n_1d_fine,   binsize_1d_fine
1564       read(1,*) n_1d_coarse, binsize_1d_coarse
1565       read(1,*) n_3d_fine,   binsize_3d_fine
1566       read(1,*) n_3d_coarse, binsize_3d_coarse
1567       read(1,*) n_3d_fine_project
1568
1569 CCC   Read momentum space track sector bins in {px,py,pz}:
1570
1571       read(1,*) n_px_bins,px_min,px_max
1572       read(1,*) n_py_bins,py_min,py_max
1573       read(1,*) n_pz_bins,pz_min,pz_max
1574
1575 CCC   Relative Chi-Square weights for track adjustment fitting:
1576
1577       read(1,*) chisq_wt_like_1d
1578       read(1,*) chisq_wt_unlike_1d
1579       read(1,*) chisq_wt_like_3d_fine
1580       read(1,*) chisq_wt_unlike_3d_fine
1581       read(1,*) chisq_wt_like_3d_coarse
1582       read(1,*) chisq_wt_unlike_3d_coarse
1583       read(1,*) chisq_wt_hist1_1
1584       read(1,*) chisq_wt_hist1_2
1585
1586       Close(unit=1)
1587       End If  !  ALICE Data I/O Option
1588
1589 C-----------------------------
1590       ELSE IF (mode.eq.2) THEN
1591 C-----------------------------
1592
1593 C     Open event generator text file, 'event_text.in,' and read it,
1594 C     noting each type of line input.  Write out a file called
1595 C     'event_line.flags' which identifies the type of information on
1596 C     each line where:
1597 C
1598 C        'EVENT:'  lines are assigned flag = 1
1599 C        'VERTEX:' lines are assigned flag = 2
1600 C        'TRACK:'  lines are assigned flag = 3
1601 C        'GENER:'  lines are assigned flag = 5
1602 C        All other lines are assigned flag = 0
1603
1604       If(ALICE .eq. 0) Then
1605       open(unit=2,type='old',access='sequential',
1606      1     name='event_text.in')
1607       open(unit=3,status='unknown',access='sequential',
1608      1     name='event_line.flags')
1609
1610 CCC   Set Event Counter:
1611
1612       event_counter = 0
1613 30    read(2,10,err=20,end=25) evg_label
1614 10    Format(A)
1615       if(evg_label .eq. event_line) then
1616          event_counter = event_counter + 1
1617          flag = 1
1618       else if(evg_label .eq. vertex_line) then
1619          flag = 2
1620       else if(evg_label .eq. track_line) then
1621          flag = 3
1622       else if(evg_label .eq. gener_line) then
1623          flag = 5
1624       else
1625          flag = 0
1626       end if
1627
1628       write(3,11) flag
1629 11    Format(1x,I1)
1630       go to 30         ! Return to S.N. 30 and read next line in file
1631 20    write(7,12) event_counter
1632 12    Format(5x,'Read error in event_text.in at event# ',I5,' - STOP')
1633       Stop
1634 25    Continue
1635       Close(unit=2)
1636       Close(unit=3)
1637       End If  !  ALICE Data I/O Option
1638
1639 C     Next, re-open the 'event_text.in' and 'event_line.flags' files
1640 C     again and read thru the entire files.  For each track, check its'
1641 C     particle ID and kinematics (pt,phi,eta) with respect to the
1642 C     selected particle ID type(s) for the run and the acceptances.
1643 C     Fill another file called, 'event_tracks.select,' which is the same
1644 C     as 'event_line.flags' except that the accepted tracks are marked
1645 C     with flag = 4.
1646 C
1647 C     NOTE:  Assume all vertices in 'event_text.in' are at microscopic
1648 C            distances (fermis) such that all particles in the event
1649 C            file are considered as primaries.  Also for each event
1650 C            the code will only accept tracks up to the limit imposed
1651 C            by trk_maxlen in the 'trk' table.
1652
1653       If(ALICE .eq. 1) Then
1654 CCC   For ALICE application do the following:
1655 CCC      Store number of events in 'n_events'
1656 CCC      Count number accepted tracks in each event, check wrt trk_maxlen
1657 CCC      Mark accepted tracks in all events
1658
1659          Call AliHbtp_GetNumberEvents(n_events)
1660          do i = 1,n_events
1661             Call AliHbtp_SetActiveEventNumber(i)
1662             track_counter = 0
1663             Call AliHbtp_GetNumberTracks(ntracks)
1664             do j = 1,ntracks
1665                Call AliHbtp_GetTrack(j,flag,px,py,pz,ge_pid)
1666                eg_pid = ge_pid
1667 CCC   Check if this track's particle ID is one to be used
1668
1669          pidok = 0
1670          accepok = 0
1671          if(pid(1).gt.0 .and. eg_pid.eq.pid(1)) pidok = 1
1672          if(pid(2).gt.0 .and. eg_pid.eq.pid(2)) pidok = 1
1673          if(pidok.eq.1 .and. eg_pid.le.part_maxlen) then
1674             mass = part_mass(eg_pid)
1675             Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
1676             if(pt.ge.pt_min .and. pt.le.pt_max .and.
1677      1         phi.ge.phi_min .and. phi.le.phi_max .and.
1678      2         eta.ge.eta_min .and. eta.le.eta_max) then
1679                if(track_counter .lt. trk_maxlen) then
1680                   accepok = 1
1681                else
1682                   write(7,62) trk_maxlen, event_counter
1683 62                Format(5x,'#tracks exceeds trk_maxlen = ',
1684      1                  I6,' for event#',I4)
1685                end if
1686             end if
1687          end if
1688
1689          if(pidok.eq.1 .and. accepok.eq.1) then
1690             track_counter = track_counter + 1
1691 C           write(*,*) '  FFF: 1 calling PutTrack j = ',j
1692             Call AliHbtp_PutTrack(j,flag4,px,py,pz,ge_pid)
1693          else
1694 C           write(*,*) '  FFF: 2 calling PutTrack j = ',j
1695             Call AliHbtp_PutTrack(j,flag0,px,py,pz,ge_pid)
1696          end if
1697          end do
1698          end do
1699
1700       Else If(ALICE .eq. 0) Then
1701
1702       open(unit=2,type='old',access='sequential',
1703      1     name='event_text.in')
1704       open(unit=3,type='old',access='sequential',
1705      1     name='event_line.flags')
1706       open(unit=4,status='unknown',access='sequential',
1707      1     name='event_tracks.select')
1708
1709 CCC   Set Event Counter:
1710
1711       event_counter = 0
1712 40    read(3,11,err=45,end=50) flag
1713       if(flag.eq.1) then
1714          event_counter = event_counter + 1
1715          track_counter = 0
1716       end if
1717
1718       if(flag.ne.3) then
1719          read(2,10) dummy
1720          write(4,11) flag
1721       else if(flag.eq.3) then
1722          read(2,41,err=46,end=50) ge_pid,px,py,pz,tid,start_v,
1723      1                            stop_v,eg_pid
1724 41    Format(7x,I6,3(1x,G12.5),4(1x,I6))
1725
1726 CCC   Check if the 'event_text.in' file includes non-zero PID
1727 CCC   values for the variable 'eg_pid'.  If this is zero, then
1728 CCC   use the ge_pid value.
1729          if(eg_pid.eq.0 .and. ge_pid.ne.0) eg_pid = ge_pid
1730
1731 CCC   Check if this track's particle ID is one to be used
1732
1733          pidok = 0
1734          accepok = 0
1735          if(pid(1).gt.0 .and. eg_pid.eq.pid(1)) pidok = 1
1736          if(pid(2).gt.0 .and. eg_pid.eq.pid(2)) pidok = 1
1737          if(pidok.eq.1 .and. eg_pid.le.part_maxlen) then
1738             mass = part_mass(eg_pid)
1739             Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
1740             if(pt.ge.pt_min .and. pt.le.pt_max .and.
1741      1         phi.ge.phi_min .and. phi.le.phi_max .and.
1742      2         eta.ge.eta_min .and. eta.le.eta_max) then
1743                if(track_counter .lt. trk_maxlen) then
1744                   accepok = 1
1745                else
1746                   write(7,621) trk_maxlen, event_counter
1747 621                Format(5x,'#tracks exceeds trk_maxlen = ',
1748      1                  I6,' for event#',I4)
1749                end if
1750             end if
1751          end if
1752
1753          if(pidok.eq.1 .and. accepok.eq.1) then
1754             track_counter = track_counter + 1
1755             write(4,11) flag4
1756          else
1757             write(4,11) flag
1758          end if
1759
1760       end if    ! End Flag=3 options
1761
1762       go to 40  ! Return to S.N. 40 and read next record
1763 45    write(7,60) event_counter
1764 60    Format(5x,'Read error in event_line.flags at event#',I5,
1765      1       ' - STOP')
1766       Stop
1767 46    write(7,61) event_counter
1768 61    Format(5x,'Read error in event_text.in (2nd pass) at event#',I5,
1769      1       ' - STOP')
1770       Stop
1771 50    Continue
1772
1773       n_events = event_counter - 1  !  Set # events in event_text.in file
1774 C                                   !  This is one less than the counter
1775 C                                   !  value since the last 'EVENT:' line is
1776 C                                   !  used to mark the End-of-File.
1777
1778       Close(unit=2)
1779       Close(unit=3)
1780       Close(unit=4)
1781
1782       End If    !  ALICE Data I/O Option
1783
1784 C-----------------------------
1785       ELSE IF(mode.eq.3) THEN
1786 C-----------------------------
1787
1788 C     Read the reference histograms for pairs, then for singles for one
1789 C     or two particle ID types.  Check switches, bins and mesh information
1790 C     to be sure the input reference histograms are compatible with the
1791 C     present run conditions.
1792
1793       open(unit=9,type='old',access='sequential',
1794      1     name='hbt_pair_reference.hist')
1795
1796       read(9,*) (check(i),i=1,3)
1797       read(9,*)  check(4),acheck(1),acheck(2)
1798       read(9,*)  check(5),acheck(3),acheck(4)
1799       read(9,*)  check(6),acheck(5),acheck(6)
1800       read(9,*) (check(i),i=7,9)
1801       read(9,*) (check(i),i=10,13)
1802       read(9,*) (acheck(i),i=7,10)
1803       read(9,*) num_pairs_like_ref, num_pairs_unlike_ref
1804
1805 CCC   Determine if the Input Reference pair histograms are compatible
1806 CCC   with the present run parameters:
1807
1808       ref_check = 1
1809       if(check(1)  .ne. n_pid_types ) ref_check = 0
1810       if(check(2)  .ne. pid(1)      ) ref_check = 0
1811       if(check(3)  .ne. pid(2)      ) ref_check = 0
1812       if(check(4)  .ne. n_pt_bins   ) ref_check = 0
1813       if(check(5)  .ne. n_phi_bins  ) ref_check = 0
1814       if(check(6)  .ne. n_eta_bins  ) ref_check = 0
1815       if(check(7)  .ne. switch_1d   ) ref_check = 0
1816       if(check(8)  .ne. switch_3d   ) ref_check = 0
1817       if(check(9)  .ne. switch_type ) ref_check = 0
1818       if(check(10) .ne. n_1d_fine   ) ref_check = 0
1819       if(check(11) .ne. n_1d_coarse ) ref_check = 0
1820       if(check(12) .ne. n_3d_fine   ) ref_check = 0
1821       if(check(13) .ne. n_3d_coarse ) ref_check = 0
1822
1823       if(abs(acheck( 1) - pt_min )         .gt. 0.000001) ref_check = 0
1824       if(abs(acheck( 2) - pt_max )         .gt. 0.000001) ref_check = 0
1825       if(abs(acheck( 3) - phi_min )        .gt. 0.000001) ref_check = 0
1826       if(abs(acheck( 4) - phi_max )        .gt. 0.000001) ref_check = 0
1827       if(abs(acheck( 5) - eta_min )        .gt. 0.000001) ref_check = 0
1828       if(abs(acheck( 6) - eta_max )        .gt. 0.000001) ref_check = 0
1829       if(abs(acheck( 7) - binsize_1d_fine) .gt. 0.000001) ref_check = 0
1830       if(abs(acheck( 8) - binsize_1d_coarse).gt.0.000001) ref_check = 0
1831       if(abs(acheck( 9) - binsize_3d_fine) .gt. 0.000001) ref_check = 0
1832       if(abs(acheck(10) - binsize_3d_coarse).gt.0.000001) ref_check = 0
1833
1834       if(ref_check .eq. 0) then
1835          write(7,70)
1836 70    Format(5x,'Reference Pair Histogram Parameters not consistent',
1837      1       ' with present run conditions - STOP')
1838          errorcode = 1
1839          Return
1840       else if(ref_check .eq. 1) then
1841
1842          if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
1843             if(switch_type.eq.1 .or. switch_type.eq.3) then
1844                read(9,*) (href_like_1d(i),i=1,n_1d_total)
1845             end if
1846             if(switch_type.eq.2 .or. switch_type.eq.3) then
1847                read(9,*) (href_unlike_1d(i),i=1,n_1d_total)
1848             end if
1849          end if   ! End 1D input option
1850
1851          if(switch_3d.gt.0) then
1852             if(switch_type.eq.1 .or. switch_type.eq.3) then
1853
1854                if(n_3d_fine.gt.0) then
1855                   do i = 1,n_3d_fine
1856                   do j = 1,n_3d_fine
1857                   do k = 1,n_3d_fine
1858                      read(9,*) href_like_3d_fine(i,j,k)
1859                   end do
1860                   end do
1861                   end do
1862                end if
1863
1864                if(n_3d_coarse.gt.0) then
1865                   do i = 1,n_3d_coarse
1866                   do j = 1,n_3d_coarse
1867                   do k = 1,n_3d_coarse
1868                      read(9,*) href_like_3d_coarse(i,j,k)
1869                   end do
1870                   end do
1871                   end do
1872                end if
1873
1874             end if
1875
1876             if(switch_type.eq.2 .or. switch_type.eq.3) then
1877
1878                if(n_3d_fine.gt.0) then
1879                   do i = 1,n_3d_fine
1880                   do j = 1,n_3d_fine
1881                   do k = 1,n_3d_fine
1882                      read(9,*) href_unlike_3d_fine(i,j,k)
1883                   end do
1884                   end do
1885                   end do
1886                end if
1887
1888                if(n_3d_coarse.gt.0) then
1889                   do i = 1,n_3d_coarse
1890                   do j = 1,n_3d_coarse
1891                   do k = 1,n_3d_coarse
1892                      read(9,*) href_unlike_3d_coarse(i,j,k)
1893                   end do
1894                   end do
1895                   end do
1896                end if
1897
1898             end if
1899
1900          end if     ! End 3D input option
1901       end if        ! End Reference Check OK/Not OK Option
1902
1903       Close(unit=9)
1904
1905 CCC   Next read the one-body histograms for 1 or 2 particle ID types:
1906
1907       open(unit=11,type='old',access='sequential',
1908      1     name='hbt_singles_reference.hist')
1909
1910       read(11,*) (check(i),i=1,3)
1911       read(11,*)  check(4),acheck(1),acheck(2)
1912       read(11,*)  check(5),acheck(3),acheck(4)
1913       read(11,*)  check(6),acheck(5),acheck(6)
1914       read(11,*) n_part_used_1_ref, n_part_used_2_ref
1915
1916 CCC   Determine if Reference one-body histograms are compatible with
1917 CCC   the present run conditions.
1918
1919       ref_check = 1
1920       if(check(1) .ne. n_pid_types)   ref_check = 0
1921       if(check(2) .ne. pid(1)     )   ref_check = 0
1922       if(check(3) .ne. pid(2)     )   ref_check = 0
1923       if(check(4) .ne. n_pt_bins  )   ref_check = 0
1924       if(check(5) .ne. n_phi_bins )   ref_check = 0
1925       if(check(6) .ne. n_eta_bins )   ref_check = 0
1926
1927       if(abs(acheck(1) - pt_min  ).gt.0.000001) ref_check = 0
1928       if(abs(acheck(2) - pt_max  ).gt.0.000001) ref_check = 0
1929       if(abs(acheck(3) - phi_min ).gt.0.000001) ref_check = 0
1930       if(abs(acheck(4) - phi_max ).gt.0.000001) ref_check = 0
1931       if(abs(acheck(5) - eta_min ).gt.0.000001) ref_check = 0
1932       if(abs(acheck(6) - eta_max ).gt.0.000001) ref_check = 0
1933
1934       if(ref_check .eq. 0) then
1935          write(7,71)
1936 71       Format(5x,'Reference One-Body Histogram parameters not ',
1937      1          'consistent with current run - STOP')
1938          errorcode = 1
1939          Return
1940       else if(ref_check .eq. 1) then
1941
1942          if(pid(1).gt.0) then
1943             read(11,*) (href1_pt_1(i) ,i=1,n_pt_bins)
1944             read(11,*) (href1_phi_1(i),i=1,n_phi_bins)
1945             read(11,*) (href1_eta_1(i),i=1,n_eta_bins)
1946          end if
1947
1948          if(pid(2).gt.0) then
1949             read(11,*) (href1_pt_2(i) ,i=1,n_pt_bins)
1950             read(11,*) (href1_phi_2(i),i=1,n_phi_bins)
1951             read(11,*) (href1_eta_2(i),i=1,n_eta_bins)
1952          end if
1953
1954       end if    ! End one-body reference histogram input
1955
1956       Close(unit=11)
1957
1958 C-----------------------------
1959       ELSE IF(mode.eq.4) THEN
1960 C-----------------------------
1961
1962 CCC   Read the next event from 'event_text.in' and load accepted tracks
1963 C     into the 'trk' data structure using the flag information about each
1964 C     line type in the file 'event_tracks.select'.
1965 C
1966 C     For this mode to run successfully the calling program must:
1967 C       (1) initially set the event_line_counter = 0
1968 C       (2) open the 'event_text.in' and 'event_tracks.select' files
1969 C           as units 2 and 4, respectively.
1970 C       (3) Close units 2 and 4 when finished.
1971
1972 CCC   Initialize accepted track counters for this new event:
1973
1974       track_counter   = 0     ! Counts all accepted tracks
1975       track_counter_1 = 0     ! Counts all accepted tracks of type pid(1)
1976       track_counter_2 = 0     ! Counts all accepted tracks of type pid(2)
1977
1978       If(ALICE .eq. 1) Then
1979          Call AliHbtp_GetNumberTracks(ntracks)
1980          do i = 1,ntracks
1981             Call AliHbtp_GetTrack(i,flag,px,py,pz,ge_pid)
1982             eg_pid = ge_pid
1983             if(flag.eq.flag4) then
1984          track_counter = track_counter + 1
1985
1986          if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
1987             track_counter_1 = track_counter_1 + 1
1988          end if
1989
1990          if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
1991             track_counter_2 = track_counter_2 + 1
1992          end if
1993
1994          mass = part_mass(eg_pid)
1995          Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
1996          trk_ge_pid(track_counter)        = eg_pid
1997          trk_px(track_counter)            = px
1998          trk_py(track_counter)            = py
1999          trk_pz(track_counter)            = pz
2000          trk_id(track_counter)            = track_counter
2001          trk_start_vertex(track_counter)  = 0
2002          trk_stop_vertex(track_counter)   = 0
2003          trk_event_line(track_counter)    = 0
2004          trk_flag(track_counter)          = 0
2005          trk_px_sec(track_counter)        = 0
2006          trk_py_sec(track_counter)        = 0
2007          trk_pz_sec(track_counter)        = 0
2008          trk_sector(track_counter)        = 0
2009          trk_out_flag(track_counter)      = 0
2010          trk_merge_flag(track_counter)    = 0
2011          trk_E(track_counter)             = E
2012          trk_pt(track_counter)            = pt
2013          trk_phi(track_counter)           = phi
2014          trk_eta(track_counter)           = eta
2015         end if
2016        end do
2017          n_part_1_trk   = track_counter_1
2018          n_part_2_trk   = track_counter_2
2019          n_part_tot_trk = track_counter
2020
2021       Else If(ALICE .eq. 0) Then
2022
2023 80    read(4,11,err=81,end=82) flag
2024       event_line_counter = event_line_counter + 1
2025
2026       if(flag .ne. 4) then
2027          read(2,10,err=83,end=82) dummy
2028       else if(flag .eq. 4) then
2029          read(2,41) ge_pid,px,py,pz,tid,start_v,stop_v,eg_pid
2030
2031 CCC   Check if the 'event_text.in' file includes non-zero PID
2032 CCC   values for the variable 'eg_pid'.  If this is zero, then
2033 CCC   use the ge_pid value.
2034          if(eg_pid.eq.0 .and. ge_pid.ne.0) eg_pid = ge_pid
2035
2036          track_counter = track_counter + 1
2037
2038          if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2039             track_counter_1 = track_counter_1 + 1
2040          end if
2041
2042          if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2043             track_counter_2 = track_counter_2 + 1
2044          end if
2045
2046          mass = part_mass(eg_pid)
2047          Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
2048          trk_ge_pid(track_counter)        = eg_pid
2049          trk_px(track_counter)            = px
2050          trk_py(track_counter)            = py
2051          trk_pz(track_counter)            = pz
2052          trk_id(track_counter)            = track_counter
2053          trk_start_vertex(track_counter)  = start_v
2054          trk_stop_vertex(track_counter)   = stop_v
2055          trk_event_line(track_counter)    = event_line_counter
2056          trk_flag(track_counter)          = 0
2057          trk_px_sec(track_counter)        = 0
2058          trk_py_sec(track_counter)        = 0
2059          trk_pz_sec(track_counter)        = 0
2060          trk_sector(track_counter)        = 0
2061          trk_out_flag(track_counter)      = 0
2062          trk_merge_flag(track_counter)    = 0
2063          trk_E(track_counter)             = E
2064          trk_pt(track_counter)            = pt
2065          trk_phi(track_counter)           = phi
2066          trk_eta(track_counter)           = eta
2067       end if
2068
2069       if(flag.ne.1) then
2070          go to 80   ! Return to S.N. 80 and read next record in file
2071       else if(flag.eq.1) then
2072          n_part_1_trk   = track_counter_1
2073          n_part_2_trk   = track_counter_2
2074          n_part_tot_trk = track_counter
2075       end if
2076
2077 82    Return
2078 81    write(7,84)
2079 84    Format(5x,'Read error from file event_tracks.select for mode=4',
2080      1       ' - STOP')
2081       Stop
2082 83    write(7,85)
2083 85    Format(5x,'Read error from file event_text.in for mode=4',
2084      1       ' - STOP')
2085       Stop
2086       End If  !  ALICE Data I/O Option
2087
2088 C-----------------------------
2089       ELSE IF(mode.eq.5) THEN
2090 C-----------------------------
2091
2092 CCC   Read the next event from 'event_text.in' and load accepted tracks
2093 C     into the 'trk2' data structure using the flag information about each
2094 C     line type in the file 'event_tracks.select'.
2095 C
2096 C     For this mode to run successfully the calling program must:
2097 C       (1) initially set the event_line_counter = 0
2098 C       (2) open the 'event_text.in' and 'event_tracks.select' files
2099 C           as units 2 and 4, respectively.
2100 C       (3) Close units 2 and 4 when finished.
2101
2102 CCC   Initialize accepted track counters for this new event:
2103
2104       track_counter   = 0     ! Counts all accepted tracks
2105       track_counter_1 = 0     ! Counts all accepted tracks of type pid(1)
2106       track_counter_2 = 0     ! Counts all accepted tracks of type pid(2)
2107
2108       If(ALICE .eq. 1) Then
2109          Call AliHbtp_GetNumberTracks(ntracks)
2110          do i = 1,ntracks
2111             Call AliHbtp_GetTrack(i,flag,px,py,pz,ge_pid)
2112             eg_pid = ge_pid
2113             if(flag.eq.flag4) then
2114          track_counter = track_counter + 1
2115
2116          if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2117             track_counter_1 = track_counter_1 + 1
2118          end if
2119
2120          if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2121             track_counter_2 = track_counter_2 + 1
2122          end if
2123
2124          mass = part_mass(eg_pid)
2125          Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
2126          trk2_ge_pid(track_counter)        = eg_pid
2127          trk2_px(track_counter)            = px
2128          trk2_py(track_counter)            = py
2129          trk2_pz(track_counter)            = pz
2130          trk2_id(track_counter)            = track_counter
2131          trk2_start_vertex(track_counter)  = 0
2132          trk2_stop_vertex(track_counter)   = 0
2133          trk2_event_line(track_counter)    = 0
2134          trk2_flag(track_counter)          = 0
2135          trk2_px_sec(track_counter)        = 0
2136          trk2_py_sec(track_counter)        = 0
2137          trk2_pz_sec(track_counter)        = 0
2138          trk2_sector(track_counter)        = 0
2139          trk2_out_flag(track_counter)      = 0
2140          trk2_merge_flag(track_counter)    = 0
2141          trk2_E(track_counter)             = E
2142          trk2_pt(track_counter)            = pt
2143          trk2_phi(track_counter)           = phi
2144          trk2_eta(track_counter)           = eta
2145         end if
2146        end do
2147          n_part_1_trk2   = track_counter_1
2148          n_part_2_trk2   = track_counter_2
2149          n_part_tot_trk2 = track_counter
2150
2151       Else If(ALICE.eq.0) Then
2152
2153 90    read(4,11,err=91,end=92) flag
2154       event_line_counter = event_line_counter + 1
2155
2156       if(flag .ne. 4) then
2157          read(2,10,err=93,end=92) dummy
2158       else if(flag .eq. 4) then
2159          read(2,41) ge_pid,px,py,pz,tid,start_v,stop_v,eg_pid
2160
2161 CCC   Check if the 'event_text.in' file includes non-zero PID
2162 CCC   values for the variable 'eg_pid'.  If this is zero, then
2163 CCC   use the ge_pid value.
2164          if(eg_pid.eq.0 .and. ge_pid.ne.0) eg_pid = ge_pid
2165
2166          track_counter = track_counter + 1
2167
2168          if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2169             track_counter_1 = track_counter_1 + 1
2170          end if
2171
2172          if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2173             track_counter_2 = track_counter_2 + 1
2174          end if
2175
2176          mass = part_mass(eg_pid)
2177          Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
2178          trk2_ge_pid(track_counter)        = eg_pid
2179          trk2_px(track_counter)            = px
2180          trk2_py(track_counter)            = py
2181          trk2_pz(track_counter)            = pz
2182          trk2_id(track_counter)            = track_counter
2183          trk2_start_vertex(track_counter)  = start_v
2184          trk2_stop_vertex(track_counter)   = stop_v
2185          trk2_event_line(track_counter)    = event_line_counter
2186          trk2_flag(track_counter)          = 0
2187          trk2_px_sec(track_counter)        = 0
2188          trk2_py_sec(track_counter)        = 0
2189          trk2_pz_sec(track_counter)        = 0
2190          trk2_sector(track_counter)        = 0
2191          trk2_out_flag(track_counter)      = 0
2192          trk2_merge_flag(track_counter)    = 0
2193          trk2_E(track_counter)             = E
2194          trk2_pt(track_counter)            = pt
2195          trk2_phi(track_counter)           = phi
2196          trk2_eta(track_counter)           = eta
2197       end if
2198
2199       if(flag.ne.1) then
2200          go to 90   ! Return to S.N. 90 and read next record in file
2201       else if(flag.eq.1) then
2202          n_part_1_trk2   = track_counter_1
2203          n_part_2_trk2   = track_counter_2
2204          n_part_tot_trk2 = track_counter
2205       end if
2206
2207 92    Return
2208 91    write(7,94)
2209 94    Format(5x,'Read error from file event_tracks.select for mode=5',
2210      1       ' - STOP')
2211       Stop
2212 93    write(7,95)
2213 95    Format(5x,'Read error from file event_text.in for mode=5',
2214      1       ' - STOP')
2215       Stop
2216
2217       End If   !  ALICE Data I/O Option
2218
2219 C-----------------------------
2220       ELSE IF(mode.eq.6) THEN
2221 C-----------------------------
2222
2223 CCC   Read finite source size Coulomb pair correlation corrections and
2224 CCC   interpolate to requested source radius and save the results for q,
2225 CCC   like and unlike pairs in common/coulomb/.
2226
2227       if(switch_coulomb.eq.3 .and. Q0.ge.coulradmin .and.
2228      1   Q0.le.coulradmax) then
2229
2230 CCC   Initially, read and interpolate like pair Coulomb corrections:
2231 C     ALICE
2232
2233       If(ALICE .eq. 1) then
2234
2235         CALL GETENVF('ALICE_ROOT',CHROOT)
2236         LNROOT = LNBLNK(CHROOT)
2237
2238         IF(LNROOT.LE.0) THEN
2239          PRINT*,'**********************************'
2240          PRINT*,'*        HBT PROCESSOR           *'
2241          PRINT*,'*        -----------             *'
2242          PRINT*,'*   DATA File  not found         *'
2243          PRINT*,'*         Program STOP           *'
2244          PRINT*,'*   Check ALICE_ROOT environment *'
2245          PRINT*,'*           variable             *'
2246          PRINT*,'**********************************'
2247          errorcode = 1
2248          return
2249         ENDIF
2250
2251         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_06.dat'
2252         open(unit=21,type='old',access='sequential',
2253      1       name=FILNAM)
2254
2255         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_08.dat'
2256         open(unit=22,type='old',access='sequential',
2257      1       name=FILNAM)
2258
2259         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_10.dat'
2260         open(unit=23,type='old',access='sequential',
2261      1       name=FILNAM)
2262
2263         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_12.dat'
2264         open(unit=24,type='old',access='sequential',
2265      1       name=FILNAM)
2266
2267         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_14.dat'
2268         open(unit=25,type='old',access='sequential',
2269      1       name=FILNAM)
2270
2271         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_16.dat'
2272         open(unit=26,type='old',access='sequential',
2273      1       name=FILNAM)
2274
2275         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_18.dat'
2276         open(unit=27,type='old',access='sequential',
2277      1       name=FILNAM)
2278
2279       ELSE
2280         open(unit=21,type='old',access='sequential',
2281      1       name='cpp_06.dat')
2282         open(unit=22,type='old',access='sequential',
2283      1       name='cpp_08.dat')
2284         open(unit=23,type='old',access='sequential',
2285      1       name='cpp_10.dat')
2286         open(unit=24,type='old',access='sequential',
2287      1       name='cpp_12.dat')
2288         open(unit=25,type='old',access='sequential',
2289      1       name='cpp_14.dat')
2290         open(unit=26,type='old',access='sequential',
2291      1       name='cpp_16.dat')
2292         open(unit=27,type='old',access='sequential',
2293      1       name='cpp_18.dat')
2294       ENDIF
2295
2296
2297       do i = 1,max_c2_coul
2298          do j = 1,ncoulradsteps
2299             read(20+j,*) q_coul(i), function(j)
2300          end do
2301          Call AliHbtp_interp(coulradmin,coulradmax,coulradstep,
2302      1        ncoulradsteps,function,20,Q0,c2_coul_like(i))
2303       end do
2304
2305       close(unit=21)
2306       close(unit=22)
2307       close(unit=23)
2308       close(unit=24)
2309       close(unit=25)
2310       close(unit=26)
2311       close(unit=27)
2312
2313 CCC   Next read and interpolate the unlike pair Coulomb corrections:
2314
2315       If(ALICE .eq. 1) then
2316         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_06.dat'
2317         open(unit=31,type='old',access='sequential',
2318      1       name=FILNAM)
2319
2320         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_08.dat'
2321         open(unit=32,type='old',access='sequential',
2322      1       name=FILNAM)
2323
2324         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_10.dat'
2325         open(unit=33,type='old',access='sequential',
2326      1       name=FILNAM)
2327
2328         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_12.dat'
2329         open(unit=34,type='old',access='sequential',
2330      1       name=FILNAM)
2331
2332         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_14.dat'
2333         open(unit=35,type='old',access='sequential',
2334      1       name=FILNAM)
2335
2336         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_16.dat'
2337         open(unit=36,type='old',access='sequential',
2338      1       name=FILNAM)
2339
2340         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_18.dat'
2341         open(unit=37,type='old',access='sequential',
2342      1       name=FILNAM)
2343
2344       else
2345         open(unit=31,type='old',access='sequential',
2346      1       name='cpm_06.dat')
2347         open(unit=32,type='old',access='sequential',
2348      1       name='cpm_08.dat')
2349         open(unit=33,type='old',access='sequential',
2350      1       name='cpm_10.dat')
2351         open(unit=34,type='old',access='sequential',
2352      1       name='cpm_12.dat')
2353         open(unit=35,type='old',access='sequential',
2354      1       name='cpm_14.dat')
2355         open(unit=36,type='old',access='sequential',
2356      1       name='cpm_16.dat')
2357         open(unit=37,type='old',access='sequential',
2358      1       name='cpm_18.dat')
2359       EndIf
2360
2361       do i = 1,max_c2_coul
2362          do j = 1,ncoulradsteps
2363             read(30+j,*) q_coul(i), function(j)
2364          end do
2365          Call AliHbtp_interp(coulradmin,coulradmax,coulradstep,
2366      1        ncoulradsteps,function,20,Q0,c2_coul_unlike(i))
2367       end do
2368
2369       close(unit=31)
2370       close(unit=32)
2371       close(unit=33)
2372       close(unit=34)
2373       close(unit=35)
2374       close(unit=36)
2375       close(unit=37)
2376
2377 CCC   Convert the input q values which are in MeV/c, to GeV/c:
2378
2379       do i = 1,max_c2_coul
2380          q_coul(i) = 0.001*q_coul(i)
2381       end do
2382
2383       end if
2384
2385
2386 C----------------------------
2387       ELSE IF(mode.eq.7) THEN
2388 C----------------------------
2389
2390 CCC   Read next event from 'event_text.in', load accepted tracks into 'trk'
2391 CCC   data structure using the flag information in the file
2392 CCC   'event_tracks.select', copy contents of 'event_text.in' and
2393 CCC   'event_tracks.select', for this one event only, into temporary files
2394 CCC   'event_text_aux.in' and 'event_tracks_aux.select', respectively.
2395 C
2396 C     For this mode to run successfully the calling program must:
2397 C       (1) initially set the event_line_counter = 0
2398 C       (2) open the 'event_text.in' and 'event_tracks.select' files
2399 C           as units 2 and 4, respectively.
2400 C       (3) Close units 2 and 4 when finished.
2401
2402 CCC   Initialize accepted track counters for this new event:
2403
2404       track_counter   = 0     ! Counts all accepted tracks
2405       track_counter_1 = 0     ! Counts all accepted tracks of type pid(1)
2406       track_counter_2 = 0     ! Counts all accepted tracks of type pid(2)
2407
2408       If(ALICE .eq. 1) Then
2409          Call AliHbtp_GetNumberTracks(ntracks)
2410          do i = 1,ntracks
2411             Call AliHbtp_GetTrack(i,flag,px,py,pz,ge_pid)
2412             eg_pid = ge_pid
2413             if(flag.eq.flag4) then
2414          track_counter = track_counter + 1
2415
2416          if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2417             track_counter_1 = track_counter_1 + 1
2418          end if
2419
2420          if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2421             track_counter_2 = track_counter_2 + 1
2422          end if
2423
2424          mass = part_mass(eg_pid)
2425          Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
2426          trk_ge_pid(track_counter)        = eg_pid
2427          trk_px(track_counter)            = px
2428          trk_py(track_counter)            = py
2429          trk_pz(track_counter)            = pz
2430          trk_id(track_counter)            = track_counter
2431          trk_start_vertex(track_counter)  = 0
2432          trk_stop_vertex(track_counter)   = 0
2433          trk_event_line(track_counter)    = 0
2434          trk_flag(track_counter)          = 0
2435          trk_px_sec(track_counter)        = 0
2436          trk_py_sec(track_counter)        = 0
2437          trk_pz_sec(track_counter)        = 0
2438          trk_sector(track_counter)        = 0
2439          trk_out_flag(track_counter)      = 0
2440          trk_merge_flag(track_counter)    = 0
2441          trk_E(track_counter)             = E
2442          trk_pt(track_counter)            = pt
2443          trk_phi(track_counter)           = phi
2444          trk_eta(track_counter)           = eta
2445         end if
2446        end do
2447          n_part_1_trk   = track_counter_1
2448          n_part_2_trk   = track_counter_2
2449          n_part_tot_trk = track_counter
2450
2451       Else If(ALICE .eq. 0) Then
2452
2453 CCC   Open temporary files:
2454
2455       open(unit=12,status='unknown',access='sequential',
2456      1     name='event_text_aux.in')
2457       open(unit=14,status='unknown',access='sequential',
2458      1     name='event_tracks_aux.select')
2459
2460 100   read(4,11,err=101,end=102) flag
2461       event_line_counter = event_line_counter + 1
2462       write(14,11) flag
2463       if(flag.eq.1) then
2464          read(2,10,err=103,end=102) comment_event_label
2465          write(12,10) comment_event_label
2466       else if(flag .eq. 2) then
2467          read(2,10,err=103,end=102) vertex_label
2468          write(12,10) vertex_label
2469       else if(flag .eq. 3) then
2470          read(2,10,err=103,end=102) comment_event_label
2471          write(12,10) comment_event_label
2472       else if(flag .eq. 5) then
2473          read(2,10,err=103,end=102) gener_label
2474          write(12,10) gener_label
2475       else if(flag .eq. 4) then
2476          read(2,41) ge_pid,px,py,pz,tid,start_v,stop_v,eg_pid
2477          write(12,41) ge_pid,px,py,pz,tid,start_v,stop_v,eg_pid
2478
2479 CCC   Check if the 'event_text.in' file includes non-zero PID
2480 CCC   values for the variable 'eg_pid'.  If this is zero, then
2481 CCC   use the ge_pid value.
2482          if(eg_pid.eq.0 .and. ge_pid.ne.0) eg_pid = ge_pid
2483
2484          track_counter = track_counter + 1
2485
2486          if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2487             track_counter_1 = track_counter_1 + 1
2488          end if
2489
2490          if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2491             track_counter_2 = track_counter_2 + 1
2492          end if
2493
2494          mass = part_mass(eg_pid)
2495          Call Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,1)
2496          trk_ge_pid(track_counter)        = eg_pid
2497          trk_px(track_counter)            = px
2498          trk_py(track_counter)            = py
2499          trk_pz(track_counter)            = pz
2500          trk_id(track_counter)            = track_counter
2501          trk_start_vertex(track_counter)  = start_v
2502          trk_stop_vertex(track_counter)   = stop_v
2503          trk_event_line(track_counter)    = event_line_counter
2504          trk_flag(track_counter)          = 0
2505          trk_px_sec(track_counter)        = 0
2506          trk_py_sec(track_counter)        = 0
2507          trk_pz_sec(track_counter)        = 0
2508          trk_sector(track_counter)        = 0
2509          trk_out_flag(track_counter)      = 0
2510          trk_merge_flag(track_counter)    = 0
2511          trk_E(track_counter)             = E
2512          trk_pt(track_counter)            = pt
2513          trk_phi(track_counter)           = phi
2514          trk_eta(track_counter)           = eta
2515       else
2516          read(2,10,err=103,end=102) comment_event_label
2517          write(12,10) comment_event_label
2518       end if
2519
2520       if(flag.ne.1) then
2521          go to 100   ! Return to S.N. 100 and read next record in file
2522       else if(flag.eq.1) then
2523          n_part_1_trk   = track_counter_1
2524          n_part_2_trk   = track_counter_2
2525          n_part_tot_trk = track_counter
2526       end if
2527
2528 102   Close(unit=12)
2529       Close(unit=14)
2530       Return
2531 101   write(7,104)
2532 104   Format(5x,'Read error from file event_tracks.select for mode=7',
2533      1       ' - STOP')
2534       Stop
2535 103   write(7,105)
2536 105   Format(5x,'Read error from file event_text.in for mode=7',
2537      1       ' - STOP')
2538       Stop
2539
2540       End If  !  ALICE Data I/O Option
2541
2542 C----------------------------
2543       ELSE IF(mode.eq.8) THEN
2544 C----------------------------
2545
2546 CCC   Read contents of 'event_text_aux.in' using the flag values in
2547 CCC   tmp. file 'event_tracks_aux.select' and copy this into the final
2548 CCC   output event file, 'event_hbt_text.out', where the momentum values
2549 CCC   of the accepted tracks in the initial input event file are replaced
2550 CCC   with the adjusted/correlated values obtained from the 'trk' table.
2551 C
2552 C     For this to work successfully the calling program must:
2553 C      (1) initially set the event_line_counter = 0
2554 C      (2) open the 'event_hbt_text.out' file as unit = 10
2555 C      (3) Close unit 10 when finished
2556
2557 CCC   Initialize accepted track counters:
2558
2559       track_counter = 0
2560
2561       If(ALICE .eq. 1) Then
2562          Call AliHbtp_GetNumberTracks(ntracks)
2563          do i = 1,ntracks
2564             Call AliHbtp_GetTrack(i,flag,px,py,pz,ge_pid)
2565             if(flag.eq.flag4) then
2566               track_counter = track_counter + 1
2567               if(trk_accep .ge. 1.000  .or. (trk_accep .lt. 1.00
2568      1           .and. hbtpran(irand) .le. trk_accep)) then
2569 C                 write(*,*) '  FFF: 3 calling PutTrack i = ',i
2570                  Call AliHbtp_PutTrack(i,flag,
2571      1                trk_px(track_counter),
2572      2                trk_py(track_counter),
2573      3                trk_pz(track_counter),
2574      4                ge_pid)
2575               end if
2576             end if
2577          end do
2578
2579       Else If(ALICE .eq. 0) Then
2580
2581 CCC   Open temporary, auxiliary files:
2582
2583       open(unit=12,type='old',access='sequential',
2584      1     name='event_text_aux.in')
2585       open(unit=14,type='old',access='sequential',
2586      1     name='event_tracks_aux.select')
2587
2588 120   read(14,11,err=121,end=122) flag
2589       file10_line_counter = file10_line_counter + 1
2590       if(flag.eq.1) then
2591          read(12,10,err=123,end=122) comment_event_label
2592          write(10,10) comment_event_label
2593       else if(flag .eq. 2) then
2594          read(12,10,err=123,end=122) vertex_label
2595          write(10,10) vertex_label
2596       else if(flag .eq. 3) then
2597          read(12,10,err=123,end=122) comment_event_label
2598          write(10,10) comment_event_label
2599       else if(flag .eq. 5) then
2600          read(12,10,err=123,end=122) gener_label
2601          write(10,10) gener_label
2602       else if(flag .eq. 4) then
2603          read(12,41,err=123,end=122)
2604      1      ge_pid,px,py,pz,tid,start_v,stop_v,eg_pid
2605          track_counter = track_counter + 1
2606          if(tid.eq.0) tid = trk_id(track_counter)
2607          if(trk_event_line(track_counter).eq.file10_line_counter)then
2608             if(trk_accep .ge. 1.000  .or. (trk_accep .lt. 1.00
2609      1         .and. hbtpran(irand) .le. trk_accep)) then
2610             write(10,841)ge_pid                   ,
2611      1                   trk_px(track_counter)    ,
2612      2                   trk_py(track_counter)    ,
2613      3                   trk_pz(track_counter)    ,
2614      4                   tid                      ,
2615      5                   start_v                  ,
2616      6                   stop_v                   ,
2617      7                   trk_ge_pid(track_counter)
2618 841         Format('TRACK:',1x,I6,3(1x,G12.5),4(1x,I6))
2619             end if
2620          else
2621             write(7,127)
2622             write(7,126) track_counter, trk_event_line(track_counter),
2623      1                file10_line_counter
2624 127         Format(5x,'Track table rows and Event file line count ',
2625      1          'out-of-synch. - STOP')
2626 126         Format(5x,'track_counter, trk().event_line,',
2627      1          'file10_line_counter = ',3I10)
2628             Stop
2629          end if
2630       else
2631          read(12,10,err=123,end=122) comment_event_label
2632          write(10,10) comment_event_label
2633       end if
2634
2635       if(flag .ne. 1) go to 120  !  Return to S.N. 120 and read next record
2636
2637 122   Close(unit=12,status='delete')
2638       Close(unit=14,status='delete')
2639       Return
2640 121   write(7,124)
2641 124   Format(5x,'Read error from file event_tracks_aux.select',
2642      1       ' for mode = 8 - STOP')
2643       Stop
2644 123   write(7,125)
2645 125   Format(5x,'Read error from file event_text_aux.in',
2646      1       ' for mode = 8 - STOP')
2647       Stop
2648
2649       End If  !  ALICE Data I/O Option
2650
2651 C-----------------
2652       END IF     ! End of read_data mode selection options
2653 C-----------------
2654
2655       Return
2656       END
2657
2658 C-----------------------------------------------------------------------
2659
2660
2661       subroutine getref_hist
2662       implicit none
2663
2664 CCC   This subroutine controls the reading or calculation and output
2665 CCC   of the several reference histograms.  These include:
2666 CCC      (a) the one-body {pt,phi,eta} 1D distributions for 1 or 2
2667 CCC          particle ID types.
2668 CCC      (b) the two-body pair-wise histograms for like and unlike
2669 CCC          pairs; for 1D and/or 3D fine mesh and 3D coarse mesh
2670 CCC          distributions.
2671
2672       Include 'common_parameters.inc'
2673       Include 'common_mesh.inc'
2674       Include 'common_histograms.inc'
2675
2676       Include 'common_track.inc'
2677       Include 'common_track2.inc'
2678       Include 'common_sec_track.inc'
2679       Include 'common_sec_track2.inc'
2680       Include 'common_particle.inc'
2681
2682 CCC   Local Variable Type Declarations:
2683
2684       integer*4 i,ipt,iphi,ieta,sign_toggle
2685
2686       if(ref_control .eq. 1) then
2687
2688 CCC   read pair and one-body reference histograms:
2689          Call read_data(3)
2690       else if(ref_control .eq. 2) then
2691
2692 CCC   calculate the pair and one-body histograms:
2693 CCC   Open event and flag files:
2694
2695       If(ALICE .eq. 0) Then
2696          open(unit=2,type='old',access='sequential',
2697      1        name='event_text.in')
2698          open(unit=4,type='old',access='sequential',
2699      1        name='event_tracks.select')
2700       End If
2701
2702 CCC   Initialize counters:
2703
2704          n_part_used_1_ref    = 0
2705          n_part_used_2_ref    = 0
2706          num_pairs_like_ref   = 0
2707          num_pairs_unlike_ref = 0
2708          event_line_counter   = 0
2709
2710 CCC   Read event header lines (no tracks are in this part):
2711       If(ALICE .eq. 0) Then
2712          Call read_data(4)
2713       End If
2714
2715 CCC   Set toggle switch to alternate between loading event tracks into
2716 CCC   table 'trk' and table 'trk2':
2717          sign_toggle = 1
2718
2719 CCC   Start Event Loop:
2720
2721          do i = 1,n_events
2722            If(ALICE .eq. 1) Then
2723               Call AliHbtp_SetActiveEventNumber(i)
2724            End If
2725            if(sign_toggle .eq. 1) then   ! Put tracks into 'trk'
2726              Call read_data(4)
2727              Call tindex(1,0)
2728              Call stm_build(1,0,0)
2729              if(pid(1) .gt. 0) then
2730                Call histog1(1,0,1,pid(1),0.0,0.0,0.0)
2731                n_part_used_1_ref = n_part_used_1_ref + n_part_used_1_trk
2732
2733                do ipt = 1,n_pt_bins
2734                href1_pt_1(ipt) = href1_pt_1(ipt) + hist1_pt_1(ipt)
2735                end do
2736
2737                do iphi = 1,n_phi_bins
2738                href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2739                end do
2740
2741                do ieta = 1,n_eta_bins
2742                href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2743                end do
2744              end if
2745
2746              if(pid(2) .gt. 0) then
2747                Call histog1(1,0,2,pid(2),0.0,0.0,0.0)
2748                n_part_used_2_ref = n_part_used_2_ref + n_part_used_2_trk
2749
2750                do ipt = 1,n_pt_bins
2751                href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2752                end do
2753
2754                do iphi = 1,n_phi_bins
2755                href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2756                end do
2757
2758                do ieta = 1,n_eta_bins
2759                href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2760                end do
2761              end if
2762
2763            else if(sign_toggle .eq. (-1)) then  ! Put tracks into 'trk2'
2764              Call read_data(5)
2765              Call tindex(2,0)
2766              Call stm_build(2,0,0)
2767              if(pid(1) .gt. 0) then
2768                Call histog1(4,0,1,pid(1),0.0,0.0,0.0)
2769                n_part_used_1_ref = n_part_used_1_ref +n_part_used_1_trk2
2770
2771                do ipt = 1,n_pt_bins
2772                href1_pt_1(ipt) = href1_pt_1(ipt) + hist1_pt_1(ipt)
2773                end do
2774
2775                do iphi = 1,n_phi_bins
2776                href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2777                end do
2778
2779                do ieta = 1,n_eta_bins
2780                href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2781                end do
2782              end if
2783
2784              if(pid(2) .gt. 0) then
2785                Call histog1(4,0,2,pid(2),0.0,0.0,0.0)
2786                n_part_used_2_ref = n_part_used_2_ref +n_part_used_2_trk2
2787
2788                do ipt = 1,n_pt_bins
2789                href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2790                end do
2791
2792                do iphi = 1,n_phi_bins
2793                href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2794                end do
2795
2796                do ieta = 1,n_eta_bins
2797                href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2798                end do
2799              end if
2800
2801            end if   !  End read and load to trk or trk2 option
2802
2803            sign_toggle = -sign_toggle
2804
2805            if(i .gt. 1) then    !  Compute 2-body reference histograms
2806               Call histog2(4,0,0,0,0,0.0,0.0,0.0,0.0)
2807               num_pairs_like_ref = num_pairs_like_ref
2808      1           + n_part_used_1_trk * n_part_used_1_trk2
2809      2           + n_part_used_2_trk * n_part_used_2_trk2
2810               num_pairs_unlike_ref = num_pairs_unlike_ref
2811      1           + n_part_used_1_trk * n_part_used_2_trk2
2812      2           + n_part_used_2_trk * n_part_used_1_trk2
2813            end if
2814
2815         end do   !  End of Event Loop
2816
2817 CCC   Write out the pair and one-body reference Histograms:
2818       Call write_data(2,0)
2819
2820       If(ALICE .eq. 0) Then
2821       Close(unit=2)
2822       Close(unit=4)
2823       End If
2824
2825       end if  !  End Reference Histogram read/calculate option
2826
2827       Return
2828       END
2829
2830 C----------------------------------------------------------------------
2831
2832
2833       subroutine AliHbtp_interp(rmin,rmax,rstep,nrsteps,function,
2834      1                          ndim,r,answer)
2835       implicit none
2836
2837 CCC   This routine interpolates the function values and puts the result
2838 CCC   into 'answer'.  It uses 2,3 or 4 mesh points which must be equally
2839 CCC   spaced.  The method uses the Lagrange interpolation formulas given
2840 CCC   in Abramowitz and Stegun, ``Handbook of Mathematical Functions,''
2841 CCC   (Dover Publications, New York, 1970); pages 878-879.
2842
2843 CCC   Definition of Variables in the Argument List:
2844
2845 CCC   rmin     = lower limit of independent variable for input function
2846 CCC   rmax     = upper limit of independent variable for input function
2847 CCC   rstep    = step size of independent variable
2848 CCC   nrsteps  = (redundant) # of input steps
2849 CCC   function(ndim) = Array of function values to be interpolated
2850 CCC   ndim     = array dimension size in calling program
2851 CCC   r        = coordinate value of independent variable where interpolation
2852 CCC              is needed.
2853 CCC   answer   = interpolated value
2854
2855 CCC   The algorithm will use the maximum number of points in the
2856 CCC   interpolation, up to a maximum of 4
2857
2858 CCC   If the requested coordinate value, r, is out-of-range, then
2859 CCC   'answer' is returned with a 0.0 value.
2860
2861 CCC   Local Variable Type Declarations:
2862
2863       integer*4 ndim, nrsteps, ik
2864
2865       real*4 rmin,rmax,rstep,r,answer,rshift,p
2866       real*4 function(ndim),w1,w2,w3,w4
2867
2868 CCC   Check Mesh:
2869
2870       if(abs(((rmax-rmin)/float(nrsteps-1))-rstep).gt.0.00001) then
2871          write(7,10) rmin,rmax,rstep,nrsteps
2872 10       Format(5x,'Interp mesh inconsistent:',3E12.5,I5,
2873      1          ' - STOP')
2874       Return
2875       end if
2876
2877 CCC   Check range:
2878
2879       if(r .lt. rmin .or. r .gt. rmax) then
2880          write(7,11) rmin,rmax,r
2881 11       Format(5x,'Interp called with r out-of-range =',3E12.5)
2882          answer = 0.0
2883          Return
2884       end if
2885
2886 CCC   Begin interpolation:
2887
2888       if(nrsteps .eq. 2) then
2889          p = (r - rmin)/rstep
2890          answer = (1.0 - p)*function(1) + p*function(2)
2891       else if(nrsteps .eq. 3) then
2892          p = (r - (rmin + rstep))/rstep
2893          answer = 0.5*p*(p-1.0)*function(1) + (1.0 - p*p)
2894      1            *function(2) + 0.5*p*(p+1.0)*function(3)
2895       else if(nrsteps .ge. 4) then
2896          rshift = r - rmin
2897
2898          if(rshift .le. rstep) then
2899             ik = 2
2900             p = (rshift - rstep)/rstep
2901          else if(rshift .ge. (rmax - rstep - rmin)) then
2902             ik = nrsteps - 2
2903             p = (rshift - (rmax - rmin - 2.0*rstep))/rstep
2904          else
2905             ik = int(rshift/rstep + 1.000001)
2906             if(ik .le. 1) ik = 2
2907             if(ik .ge. (nrsteps-1)) ik = nrsteps - 2
2908             p = (rshift - float(ik-1)*rstep)/rstep
2909          end if
2910
2911          w1 = -p*(p-1.0)*(p-2.0)/6.0
2912          w2 = (p*p-1.0)*(p-2.0)/2.0
2913          w3 = -p*(p+1.0)*(p-2.0)/2.0
2914          w4 = p*(p*p-1.0)/6.0
2915
2916          answer = w1*function(ik-1) + w2*function(ik)
2917      1          + w3*function(ik+1) + w4*function(ik+2)
2918       end if  !  End # interplation points option
2919
2920       Return
2921       END
2922
2923 C--------------------------------------------------------------------
2924
2925
2926       subroutine Hbtp_particle_prop
2927       implicit none
2928
2929 CCC   Fill particle properties table /particle/ with Geant 3 particle ID
2930 CCC   numbers, charge (in units of |e|), mass in GeV/c**2 and lifetime
2931 CCC   in seconds.  See the Geant 3.15 Manual User's Guide, pages: CONS
2932 CCC   300-1 and -2.
2933
2934       Include 'common_particle.inc'
2935
2936 CCC   Local Variable Type Declarations:
2937
2938       integer*4 i
2939
2940       do i = 1,part_maxlen
2941       part_id(i) = i
2942       end do
2943
2944 CCC   Set Particle Masses:
2945
2946       part_mass( 1) = 0.0                   ! Gamma
2947       part_mass( 2) = 0.00051099906         ! Positron
2948       part_mass( 3) = 0.00051099906         ! Electron
2949       part_mass( 4) = 0.0                   ! Neutrino
2950       part_mass( 5) = 0.105658389           ! Muon+
2951       part_mass( 6) = 0.105658389           ! Muon-
2952       part_mass( 7) = 0.1349743             ! Pion0
2953       part_mass( 8) = 0.1395679             ! Pion+
2954       part_mass( 9) = 0.1395679             ! Pion-
2955       part_mass(10) = 0.497671              ! Kaon 0 long
2956       part_mass(11) = 0.493646              ! Kaon+
2957       part_mass(12) = 0.493646              ! Kaon-
2958       part_mass(13) = 0.93956563            ! Neutron
2959       part_mass(14) = 0.93827231            ! Proton
2960       part_mass(15) = 0.93827231            ! Antiproton
2961       part_mass(16) = 0.497671              ! Kaon 0 short
2962       part_mass(17) = 0.54745               ! Eta
2963       part_mass(18) = 1.11563               ! Lambda
2964       part_mass(19) = 1.18937               ! Sigma+
2965       part_mass(20) = 1.19255               ! Sigma0
2966       part_mass(21) = 1.197465              ! Sigma-
2967       part_mass(22) = 1.31485               ! Xi 0
2968       part_mass(23) = 1.32133               ! Xi -
2969       part_mass(24) = 1.67243               ! Omega
2970       part_mass(25) = 0.93956563            ! Antineutron
2971       part_mass(26) = 1.11563               ! Antilambda
2972       part_mass(27) = 1.18937               ! Anti-Sigma -
2973       part_mass(28) = 1.19255               ! Anti-Sigma 0
2974       part_mass(29) = 1.197465              ! Anti-Sigma +
2975       part_mass(30) = 1.31485               ! AntiXi 0
2976       part_mass(31) = 1.32133               ! AntiXi +
2977       part_mass(32) = 1.67243               ! Anti-Omega +
2978       part_mass(33) = 0.0
2979       part_mass(34) = 0.0
2980       part_mass(35) = 0.0
2981       part_mass(36) = 0.0
2982       part_mass(37) = 0.0
2983       part_mass(38) = 0.0
2984       part_mass(39) = 0.0
2985       part_mass(40) = 0.0
2986       part_mass(41) = 0.0
2987       part_mass(42) = 0.0
2988       part_mass(43) = 0.0
2989       part_mass(44) = 0.0
2990       part_mass(45) = 1.875613              ! Deuteron
2991       part_mass(46) = 2.80925               ! Triton
2992       part_mass(47) = 3.727417              ! Alpha
2993       part_mass(48) = 0.0                   ! Geantino (Fake particle)
2994       part_mass(49) = 2.80923               ! He3
2995       part_mass(50) = 0.0                   ! Cerenkov (Fake particle)
2996
2997 CCC   Set Particle Charge:
2998
2999       part_charge( 1) =  0      ! Gamma
3000       part_charge( 2) =  1      ! Positron
3001       part_charge( 3) = -1      ! Electron
3002       part_charge( 4) =  0      ! Neutrino
3003       part_charge( 5) =  1      ! Muon+
3004       part_charge( 6) = -1      ! Muon-
3005       part_charge( 7) =  0      ! Pion0
3006       part_charge( 8) =  1      ! Pion+
3007       part_charge( 9) = -1      ! Pion-
3008       part_charge(10) =  0      ! Kaon 0 long
3009       part_charge(11) =  1      ! Kaon+
3010       part_charge(12) = -1      ! Kaon-
3011       part_charge(13) =  0      ! Neutron
3012       part_charge(14) =  1      ! Proton
3013       part_charge(15) = -1      ! Antiproton
3014       part_charge(16) =  0      ! Kaon 0 short
3015       part_charge(17) =  0      ! Eta
3016       part_charge(18) =  0      ! Lambda
3017       part_charge(19) =  1      ! Sigma+
3018       part_charge(20) =  0      ! Sigma0
3019       part_charge(21) = -1      ! Sigma-
3020       part_charge(22) =  0      ! Xi 0
3021       part_charge(23) = -1      ! Xi -
3022       part_charge(24) = -1      ! Omega
3023       part_charge(25) =  0      ! Antineutron
3024       part_charge(26) =  0      ! Antilambda
3025       part_charge(27) = -1      ! Anti-Sigma -
3026       part_charge(28) =  0      ! Anti-Sigma 0
3027       part_charge(29) =  1      ! Anti-Sigma +
3028       part_charge(30) =  0      ! AntiXi 0
3029       part_charge(31) =  1      ! AntiXi +
3030       part_charge(32) =  1      ! Anti-Omega +
3031       part_charge(33) =  0
3032       part_charge(34) =  0
3033       part_charge(35) =  0
3034       part_charge(36) =  0
3035       part_charge(37) =  0
3036       part_charge(38) =  0
3037       part_charge(39) =  0
3038       part_charge(40) =  0
3039       part_charge(41) =  0
3040       part_charge(42) =  0
3041       part_charge(43) =  0
3042       part_charge(44) =  0
3043       part_charge(45) =  1      ! Deuteron
3044       part_charge(46) =  1      ! Triton
3045       part_charge(47) =  2      ! Alpha
3046       part_charge(48) =  0      ! Geantino (Fake particle)
3047       part_charge(49) =  2      ! He3
3048       part_charge(50) =  0      ! Cerenkov (Fake particle)
3049
3050 CCC   Set Particle Lifetimes:
3051
3052       part_lifetime( 1) = 1.0E+30       ! Gamma
3053       part_lifetime( 2) = 1.0E+30       ! Positron
3054       part_lifetime( 3) = 1.0E+30       ! Electron
3055       part_lifetime( 4) = 1.0E+30       ! Neutrino
3056       part_lifetime( 5) = 2.19703E-06   ! Muon+
3057       part_lifetime( 6) = 2.19703E-06   ! Muon-
3058       part_lifetime( 7) = 8.4E-17       ! Pion0
3059       part_lifetime( 8) = 2.603E-08     ! Pion+
3060       part_lifetime( 9) = 2.603E-08     ! Pion-
3061       part_lifetime(10) = 5.16E-08      ! Kaon 0 long
3062       part_lifetime(11) = 1.237E-08     ! Kaon+
3063       part_lifetime(12) = 1.237E-08     ! Kaon-
3064       part_lifetime(13) = 889.1         ! Neutron
3065       part_lifetime(14) = 1.0E+30       ! Proton
3066       part_lifetime(15) = 1.0E+30       ! Antiproton
3067       part_lifetime(16) = 8.922E-11     ! Kaon 0 short
3068       part_lifetime(17) = 5.53085E-19   ! Eta
3069       part_lifetime(18) = 2.632E-10     ! Lambda
3070       part_lifetime(19) = 7.99E-11      ! Sigma+
3071       part_lifetime(20) = 7.40E-20      ! Sigma0
3072       part_lifetime(21) = 1.479E-10     ! Sigma-
3073       part_lifetime(22) = 2.90E-10      ! Xi 0
3074       part_lifetime(23) = 1.639E-10     ! Xi -
3075       part_lifetime(24) = 8.22E-11      ! Omega
3076       part_lifetime(25) = 889.1         ! Antineutron
3077       part_lifetime(26) = 2.632E-10     ! Antilambda
3078       part_lifetime(27) = 7.99E-11      ! Anti-Sigma -
3079       part_lifetime(28) = 7.40E-20      ! Anti-Sigma 0
3080       part_lifetime(29) = 1.479E-10     ! Anti-Sigma +
3081       part_lifetime(30) = 2.90E-10      ! AntiXi 0
3082       part_lifetime(31) = 1.639E-10     ! AntiXi +
3083       part_lifetime(32) = 8.22E-11      ! Anti-Omega +
3084       part_lifetime(33) = 0.0
3085       part_lifetime(34) = 0.0
3086       part_lifetime(35) = 0.0
3087       part_lifetime(36) = 0.0
3088       part_lifetime(37) = 0.0
3089       part_lifetime(38) = 0.0
3090       part_lifetime(39) = 0.0
3091       part_lifetime(40) = 0.0
3092       part_lifetime(41) = 0.0
3093       part_lifetime(42) = 0.0
3094       part_lifetime(43) = 0.0
3095       part_lifetime(44) = 0.0
3096       part_lifetime(45) = 1.0E+30       ! Deuteron
3097       part_lifetime(46) = 1.0E+30       ! Triton
3098       part_lifetime(47) = 1.0E+30       ! Alpha
3099       part_lifetime(48) = 1.0E+30       ! Geantino (Fake particle)
3100       part_lifetime(49) = 1.0E+30       ! He3
3101       part_lifetime(50) = 1.0E+30       ! Cerenkov (Fake particle)
3102
3103       Return
3104       END
3105
3106 C----------------------------------------------------------------
3107
3108
3109       subroutine correl_model
3110       implicit none
3111
3112 CCC   This subroutine computes the requested 2-body model correlation
3113 CCC   function which is to be fitted by the track adjustment procedure.
3114 CCC   The model values are calculated on the requested fine and coarse
3115 CCC   mesh grid in momentum space.  The model values are computed at the
3116 CCC   mid point of each 1D bin or at the center of each 3D cell.  This
3117 CCC   could be refined at a later date to correspond to the integral of
3118 CCC   the model function over the bin width (cell volume) divided by the
3119 CCC   the bin width (cell volume).
3120
3121 C     The model includes the following options which are selected by the
3122 C     'switch*' parameters in common/parameters/:
3123 C
3124 C        switch_1d:          1D model as function of either Qinvar, Qtotal
3125 C                            or Q-vector
3126 C        switch_3d:          3D model as function of either the Bertsch-Pratt
3127 C                            side-out-long kinematics (but no cross term) or
3128 C                            the Yano-Koonin-Podgoretski perp-parallel-time
3129 C                            kinematics.
3130 C        switch_type:        Like and/or Unlike particles
3131 C        switch_coherence:   Purely incoherent source or a mixed incoherent-
3132 C                            coherent source.
3133 C        switch_coulomb:     Either (a) no Coulomb correction, (b) Gamow
3134 C                            factor, (c) NA35 parametrization, or (d) Pratt
3135 C                            Coulomb wave function integration for finite
3136 C                            size, spherical source.
3137 C        switch_fermi_bose:  Fermion or boson identical pairs.
3138
3139       Include 'common_parameters.inc'
3140       Include 'common_mesh.inc'
3141       Include 'common_correlations.inc'
3142
3143 CCC   Local Variable Type Declarations:
3144
3145       integer*4 i,j,k
3146
3147       real*4 R_1dsq, Rsidesq, Routsq, Rlongsq
3148       real*4 Rperpsq, Rparallelsq, R0sq
3149       real*4 sqrtlambda,fermi_bose_sign,coulomb_factor,coherence_fac
3150       real*4 q,q1,q2,q3
3151       real*4 b,b1,b2,b3
3152       real*4 massavg
3153
3154 CCC   Set Constants:
3155
3156       sqrtlambda  = sqrt(abs(lambda))
3157       R_1dsq      = ((R_1d     /hbc)**2)/2.0
3158       Rsidesq     = ((Rside    /hbc)**2)/2.0
3159       Routsq      = ((Rout     /hbc)**2)/2.0
3160       Rlongsq     = ((Rlong    /hbc)**2)/2.0
3161       Rperpsq     = ((Rperp    /hbc)**2)/2.0
3162       Rparallelsq = ((Rparallel/hbc)**2)/2.0
3163       R0sq        = ((R0       /hbc)**2)/2.0
3164
3165       fermi_bose_sign = float(switch_fermi_bose)
3166       coherence_fac   = switch_coherence*2.0*sqrtlambda*
3167      1                  (1.0 - sqrtlambda)
3168
3169 CCC   Determine average particle pair mass for Coulomb correction:
3170
3171       massavg = 0.14
3172       if(mass1.eq.0.0 .and. mass2.gt.0.0) massavg = mass2
3173       if(mass1.gt.0.0 .and. mass2.eq.0.0) massavg = mass1
3174       if(mass1.gt.0.0 .and. mass2.gt.0.0) massavg = 0.5*(mass1+mass2)
3175
3176 CCC   Compute 1D correlation model arrays:
3177
3178       If(switch_1d .ge. 1) then
3179         If(n_1d_fine .gt. 0) then        ! Fill the 1D fine mesh bins
3180           q = -0.5*binsize_1d_fine
3181           do i = 1,n_1d_fine
3182             q = q + binsize_1d_fine
3183             b = exp(-q*q*R_1dsq)
3184
3185             if(switch_type.eq.1 .or. switch_type.eq.3) then
3186               c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3187      1                         *b*b + coherence_fac*b)
3188               if(switch_coulomb.eq.0) then
3189                 coulomb_factor = 1.0
3190               else if(switch_coulomb.gt.0) then
3191                 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3192      1                       coulomb_factor)
3193               end if
3194               c2mod_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3195             end if
3196
3197             if(switch_type.eq.2 .or. switch_type.eq.3) then
3198              c2mod_unlike_1d(i) = 1.0
3199              if(switch_coulomb.eq.0) then
3200                coulomb_factor = 1.0
3201              else if(switch_coulomb.gt.0) then
3202                Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3203      1                      coulomb_factor)
3204              end if
3205              c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3206             end if
3207
3208           end do    ! End of 1D fine mesh filling do-loop
3209         end if      ! End of 1D fine mesh option
3210
3211         If(n_1d_coarse .gt. 0) then      ! Fill the 1D coarse mesh bins
3212           q = qmid_1d -0.5*binsize_1d_coarse
3213           do i = n_1d_fine + 1, n_1d_total
3214             q = q + binsize_1d_coarse
3215             b = exp(-q*q*R_1dsq)
3216
3217             if(switch_type.eq.1 .or. switch_type.eq.3) then
3218               c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3219      1                         *b*b + coherence_fac*b)
3220               if(switch_coulomb.eq.0) then
3221                 coulomb_factor = 1.0
3222               else if(switch_coulomb.gt.0) then
3223                 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3224      1                       coulomb_factor)
3225               end if
3226               c2mod_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3227             end if
3228
3229             if(switch_type.eq.2 .or. switch_type.eq.3) then
3230              c2mod_unlike_1d(i) = 1.0
3231              if(switch_coulomb.eq.0) then
3232                coulomb_factor = 1.0
3233              else if(switch_coulomb.gt.0) then
3234                Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3235      1                      coulomb_factor)
3236              end if
3237              c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3238             end if
3239
3240           end do    ! End of 1D coarse mesh filling do-loop
3241         end if      ! End of 1D coarse mesh option
3242       end if        ! End of 1D option
3243
3244 CCC   Compute 3D correlation model arrays:
3245
3246       If(switch_3d .ge. 1) Then
3247        If(n_3d_fine .gt. 0) then     !  Fill the 3D fine mesh bins
3248         q1 = -0.5*binsize_3d_fine
3249         do i = 1,n_3d_fine
3250          q1 = q1 + binsize_3d_fine
3251          if(switch_3d.eq.1) b1=exp(-q1*q1*Rsidesq)
3252          if(switch_3d.eq.2) b1=exp(-q1*q1*Rperpsq)
3253
3254           q2 = -0.5*binsize_3d_fine
3255           do j = 1,n_3d_fine
3256            q2 = q2 + binsize_3d_fine
3257            if(switch_3d.eq.1) b2=exp(-q2*q2*Routsq)
3258            if(switch_3d.eq.2) b2=exp(-q2*q2*Rparallelsq)
3259
3260             q3 = -0.5*binsize_3d_fine
3261             do k = 1,n_3d_fine
3262              q3 = q3 + binsize_3d_fine
3263              if(switch_3d.eq.1) b3=exp(-q3*q3*Rlongsq)
3264              if(switch_3d.eq.2) b3=exp(-q3*q3*R0sq)
3265
3266              b = b1*b2*b3
3267              if(switch_3d.eq.1) q = sqrt(q1*q1+q2*q2+q3*q3)
3268              if(switch_3d.eq.2) q = sqrt(q1*q1+q2*q2)