Moving required CMake version from 2.8.4 to 2.8.8
[u/mrichter/AliRoot.git] / THbtp / 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
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
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     file='hbt_log.out')
622       open(unit=8,status='unknown',access='sequential',
623      1     file='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,status='old',access='sequential',
1522      1        file='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,status='old',access='sequential',
1606      1     file='event_text.in')
1607       open(unit=3,status='unknown',access='sequential',
1608      1     file='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,status='old',access='sequential',
1703      1     file='event_text.in')
1704       open(unit=3,status='old',access='sequential',
1705      1     file='event_line.flags')
1706       open(unit=4,status='unknown',access='sequential',
1707      1     file='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,status='old',access='sequential',
1794      1     file='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,status='old',access='sequential',
1908      1     file='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,status='old',access='sequential',
2253      1       file=FILNAM)
2254
2255         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_08.dat' 
2256         open(unit=22,status='old',access='sequential',
2257      1       file=FILNAM)
2258
2259         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_10.dat' 
2260         open(unit=23,status='old',access='sequential',
2261      1       file=FILNAM)
2262
2263         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_12.dat' 
2264         open(unit=24,status='old',access='sequential',
2265      1       file=FILNAM)
2266
2267         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_14.dat' 
2268         open(unit=25,status='old',access='sequential',
2269      1       file=FILNAM)
2270
2271         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_16.dat' 
2272         open(unit=26,status='old',access='sequential',
2273      1       file=FILNAM)
2274
2275         FILNAM=CHROOT(1:LNROOT)//'/data/cpp_18.dat' 
2276         open(unit=27,status='old',access='sequential',
2277      1       file=FILNAM)
2278         
2279       ELSE
2280         open(unit=21,status='old',access='sequential',
2281      1       file='cpp_06.dat')
2282         open(unit=22,status='old',access='sequential',
2283      1       file='cpp_08.dat')
2284         open(unit=23,status='old',access='sequential',
2285      1       file='cpp_10.dat')
2286         open(unit=24,status='old',access='sequential',
2287      1       file='cpp_12.dat')
2288         open(unit=25,status='old',access='sequential',
2289      1       file='cpp_14.dat')
2290         open(unit=26,status='old',access='sequential',
2291      1       file='cpp_16.dat')
2292         open(unit=27,status='old',access='sequential',
2293      1       file='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,status='old',access='sequential',
2318      1       file=FILNAM)
2319
2320         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_08.dat' 
2321         open(unit=32,status='old',access='sequential',
2322      1       file=FILNAM)
2323
2324         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_10.dat' 
2325         open(unit=33,status='old',access='sequential',
2326      1       file=FILNAM)
2327
2328         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_12.dat' 
2329         open(unit=34,status='old',access='sequential',
2330      1       file=FILNAM)
2331
2332         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_14.dat' 
2333         open(unit=35,status='old',access='sequential',
2334      1       file=FILNAM)
2335
2336         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_16.dat' 
2337         open(unit=36,status='old',access='sequential',
2338      1       file=FILNAM)
2339
2340         FILNAM=CHROOT(1:LNROOT)//'/data/cpm_18.dat' 
2341         open(unit=37,status='old',access='sequential',
2342      1       file=FILNAM)
2343     
2344       else
2345         open(unit=31,status='old',access='sequential',
2346      1       file='cpm_06.dat')
2347         open(unit=32,status='old',access='sequential',
2348      1       file='cpm_08.dat')
2349         open(unit=33,status='old',access='sequential',
2350      1       file='cpm_10.dat')
2351         open(unit=34,status='old',access='sequential',
2352      1       file='cpm_12.dat')
2353         open(unit=35,status='old',access='sequential',
2354      1       file='cpm_14.dat')
2355         open(unit=36,status='old',access='sequential',
2356      1       file='cpm_16.dat')
2357         open(unit=37,status='old',access='sequential',
2358      1       file='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     file='event_text_aux.in')
2457       open(unit=14,status='unknown',access='sequential',
2458      1     file='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,status='old',access='sequential',
2584      1     file='event_text_aux.in')
2585       open(unit=14,status='old',access='sequential',
2586      1     file='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,status='old',access='sequential',
2697      1        file='event_text.in')
2698          open(unit=4,status='old',access='sequential',
2699      1        file='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 C         write(*,*) 'REF HISTO N Ev = ', n_events
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                if (href1_pt_1(ipt) .lt. 0 ) then
2736                  write(*,*) 'href1_pt_1 bin ',ipt,'is less then 0'
2737                endif 
2738                end do
2739
2740                do iphi = 1,n_phi_bins
2741                href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2742                if (href1_phi_1(iphi) .lt. 0 ) then
2743                  write(*,*) 'href1_phi_1 bin ',iphi,'is less then 0'
2744                endif 
2745                end do
2746
2747                do ieta = 1,n_eta_bins
2748                href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2749                if (href1_eta_1(ieta) .lt. 0 ) then
2750                  write(*,*) 'href1_eta_1 bin ',ieta,'is less then 0'
2751                endif 
2752                end do
2753              end if
2754
2755              if(pid(2) .gt. 0) then
2756                Call histog1(1,0,2,pid(2),0.0,0.0,0.0)
2757                n_part_used_2_ref = n_part_used_2_ref + n_part_used_2_trk
2758
2759                do ipt = 1,n_pt_bins
2760                href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2761                end do
2762
2763                do iphi = 1,n_phi_bins
2764                href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2765                end do
2766
2767                do ieta = 1,n_eta_bins
2768                href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2769                end do
2770              end if
2771
2772            else if(sign_toggle .eq. (-1)) then  ! Put tracks into 'trk2'
2773              Call read_data(5)
2774              Call tindex(2,0)
2775              Call stm_build(2,0,0)
2776              if(pid(1) .gt. 0) then
2777                Call histog1(4,0,1,pid(1),0.0,0.0,0.0)
2778                n_part_used_1_ref = n_part_used_1_ref +n_part_used_1_trk2
2779
2780                do ipt = 1,n_pt_bins
2781                href1_pt_1(ipt) = href1_pt_1(ipt) + hist1_pt_1(ipt)
2782                end do
2783
2784                do iphi = 1,n_phi_bins
2785                href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2786                end do
2787
2788                do ieta = 1,n_eta_bins
2789                href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2790                end do
2791              end if
2792
2793              if(pid(2) .gt. 0) then
2794                Call histog1(4,0,2,pid(2),0.0,0.0,0.0)
2795                n_part_used_2_ref = n_part_used_2_ref +n_part_used_2_trk2
2796
2797                do ipt = 1,n_pt_bins
2798                href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2799                end do
2800
2801                do iphi = 1,n_phi_bins
2802                href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2803                end do
2804
2805                do ieta = 1,n_eta_bins
2806                href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2807                end do
2808              end if
2809
2810            end if   !  End read and load to trk or trk2 option
2811
2812            sign_toggle = -sign_toggle
2813
2814            if(i .gt. 1) then    !  Compute 2-body reference histograms
2815               Call histog2(4,0,0,0,0,0.0,0.0,0.0,0.0)
2816               num_pairs_like_ref = num_pairs_like_ref
2817      1           + n_part_used_1_trk * n_part_used_1_trk2
2818      2           + n_part_used_2_trk * n_part_used_2_trk2
2819               num_pairs_unlike_ref = num_pairs_unlike_ref
2820      1           + n_part_used_1_trk * n_part_used_2_trk2
2821      2           + n_part_used_2_trk * n_part_used_1_trk2
2822 C              write(*,*) 'num_pairs_like_ref',num_pairs_like_ref
2823 C              write(*,*) 'num_pairs_unlike_ref',num_pairs_unlike_ref
2824            end if
2825
2826         end do   !  End of Event Loop
2827
2828 CCC   Write out the pair and one-body reference Histograms:
2829       Call write_data(2,0)
2830
2831       If(ALICE .eq. 0) Then
2832       Close(unit=2)
2833       Close(unit=4)
2834       End If
2835
2836       end if  !  End Reference Histogram read/calculate option
2837
2838       Return
2839       END
2840
2841 C----------------------------------------------------------------------
2842
2843
2844       subroutine AliHbtp_interp(rmin,rmax,rstep,nrsteps,function,
2845      1                          ndim,r,answer)
2846       implicit none
2847
2848 CCC   This routine interpolates the function values and puts the result
2849 CCC   into 'answer'.  It uses 2,3 or 4 mesh points which must be equally
2850 CCC   spaced.  The method uses the Lagrange interpolation formulas given
2851 CCC   in Abramowitz and Stegun, ``Handbook of Mathematical Functions,''
2852 CCC   (Dover Publications, New York, 1970); pages 878-879.
2853
2854 CCC   Definition of Variables in the Argument List:
2855
2856 CCC   rmin     = lower limit of independent variable for input function
2857 CCC   rmax     = upper limit of independent variable for input function
2858 CCC   rstep    = step size of independent variable
2859 CCC   nrsteps  = (redundant) # of input steps
2860 CCC   function(ndim) = Array of function values to be interpolated
2861 CCC   ndim     = array dimension size in calling program
2862 CCC   r        = coordinate value of independent variable where interpolation
2863 CCC              is needed.
2864 CCC   answer   = interpolated value
2865
2866 CCC   The algorithm will use the maximum number of points in the
2867 CCC   interpolation, up to a maximum of 4
2868
2869 CCC   If the requested coordinate value, r, is out-of-range, then
2870 CCC   'answer' is returned with a 0.0 value.
2871
2872 CCC   Local Variable Type Declarations:
2873
2874       integer*4 ndim, nrsteps, ik
2875
2876       real*4 rmin,rmax,rstep,r,answer,rshift,p
2877       real*4 function(ndim),w1,w2,w3,w4
2878
2879 CCC   Check Mesh:
2880
2881       if(abs(((rmax-rmin)/float(nrsteps-1))-rstep).gt.0.00001) then
2882          write(7,10) rmin,rmax,rstep,nrsteps
2883 10       Format(5x,'Interp mesh inconsistent:',3E12.5,I5,
2884      1          ' - STOP')
2885       Return
2886       end if
2887
2888 CCC   Check range:
2889      
2890       if(r .lt. rmin .or. r .gt. rmax) then
2891          write(7,11) rmin,rmax,r
2892 11       Format(5x,'Interp called with r out-of-range =',3E12.5)
2893          answer = 0.0
2894          Return
2895       end if
2896
2897 CCC   Begin interpolation:
2898
2899       if(nrsteps .eq. 2) then
2900          p = (r - rmin)/rstep
2901          answer = (1.0 - p)*function(1) + p*function(2)
2902       else if(nrsteps .eq. 3) then
2903          p = (r - (rmin + rstep))/rstep
2904          answer = 0.5*p*(p-1.0)*function(1) + (1.0 - p*p)
2905      1            *function(2) + 0.5*p*(p+1.0)*function(3)
2906       else if(nrsteps .ge. 4) then
2907          rshift = r - rmin
2908
2909          if(rshift .le. rstep) then
2910             ik = 2
2911             p = (rshift - rstep)/rstep
2912          else if(rshift .ge. (rmax - rstep - rmin)) then
2913             ik = nrsteps - 2
2914             p = (rshift - (rmax - rmin - 2.0*rstep))/rstep
2915          else
2916             ik = int(rshift/rstep + 1.000001)
2917             if(ik .le. 1) ik = 2
2918             if(ik .ge. (nrsteps-1)) ik = nrsteps - 2
2919             p = (rshift - float(ik-1)*rstep)/rstep
2920          end if
2921
2922          w1 = -p*(p-1.0)*(p-2.0)/6.0
2923          w2 = (p*p-1.0)*(p-2.0)/2.0
2924          w3 = -p*(p+1.0)*(p-2.0)/2.0
2925          w4 = p*(p*p-1.0)/6.0
2926
2927          answer = w1*function(ik-1) + w2*function(ik)
2928      1          + w3*function(ik+1) + w4*function(ik+2)
2929       end if  !  End # interplation points option
2930
2931       Return
2932       END
2933
2934 C--------------------------------------------------------------------
2935
2936
2937       subroutine Hbtp_particle_prop
2938       implicit none
2939
2940 CCC   Fill particle properties table /particle/ with Geant 3 particle ID
2941 CCC   numbers, charge (in units of |e|), mass in GeV/c**2 and lifetime
2942 CCC   in seconds.  See the Geant 3.15 Manual User's Guide, pages: CONS
2943 CCC   300-1 and -2.
2944
2945       Include 'common_particle.inc'
2946
2947 CCC   Local Variable Type Declarations:
2948
2949       integer*4 i
2950
2951       do i = 1,part_maxlen
2952       part_id(i) = i
2953       end do
2954
2955 CCC   Set Particle Masses:
2956
2957       part_mass( 1) = 0.0                   ! Gamma
2958       part_mass( 2) = 0.00051099906         ! Positron
2959       part_mass( 3) = 0.00051099906         ! Electron
2960       part_mass( 4) = 0.0                   ! Neutrino
2961       part_mass( 5) = 0.105658389           ! Muon+
2962       part_mass( 6) = 0.105658389           ! Muon-
2963       part_mass( 7) = 0.1349743             ! Pion0
2964       part_mass( 8) = 0.1395679             ! Pion+
2965       part_mass( 9) = 0.1395679             ! Pion-
2966       part_mass(10) = 0.497671              ! Kaon 0 long
2967       part_mass(11) = 0.493646              ! Kaon+
2968       part_mass(12) = 0.493646              ! Kaon-
2969       part_mass(13) = 0.93956563            ! Neutron
2970       part_mass(14) = 0.93827231            ! Proton
2971       part_mass(15) = 0.93827231            ! Antiproton
2972       part_mass(16) = 0.497671              ! Kaon 0 short
2973       part_mass(17) = 0.54745               ! Eta
2974       part_mass(18) = 1.11563               ! Lambda
2975       part_mass(19) = 1.18937               ! Sigma+
2976       part_mass(20) = 1.19255               ! Sigma0
2977       part_mass(21) = 1.197465              ! Sigma-
2978       part_mass(22) = 1.31485               ! Xi 0
2979       part_mass(23) = 1.32133               ! Xi -
2980       part_mass(24) = 1.67243               ! Omega
2981       part_mass(25) = 0.93956563            ! Antineutron
2982       part_mass(26) = 1.11563               ! Antilambda
2983       part_mass(27) = 1.18937               ! Anti-Sigma -
2984       part_mass(28) = 1.19255               ! Anti-Sigma 0
2985       part_mass(29) = 1.197465              ! Anti-Sigma +
2986       part_mass(30) = 1.31485               ! AntiXi 0
2987       part_mass(31) = 1.32133               ! AntiXi +
2988       part_mass(32) = 1.67243               ! Anti-Omega +
2989       part_mass(33) = 0.0
2990       part_mass(34) = 0.0
2991       part_mass(35) = 0.0
2992       part_mass(36) = 0.0
2993       part_mass(37) = 0.0
2994       part_mass(38) = 0.0
2995       part_mass(39) = 0.0
2996       part_mass(40) = 0.0
2997       part_mass(41) = 0.0
2998       part_mass(42) = 0.0
2999       part_mass(43) = 0.0
3000       part_mass(44) = 0.0
3001       part_mass(45) = 1.875613              ! Deuteron
3002       part_mass(46) = 2.80925               ! Triton
3003       part_mass(47) = 3.727417              ! Alpha
3004       part_mass(48) = 0.0                   ! Geantino (Fake particle)
3005       part_mass(49) = 2.80923               ! He3
3006       part_mass(50) = 0.0                   ! Cerenkov (Fake particle)
3007
3008 CCC   Set Particle Charge:
3009
3010       part_charge( 1) =  0      ! Gamma
3011       part_charge( 2) =  1      ! Positron
3012       part_charge( 3) = -1      ! Electron
3013       part_charge( 4) =  0      ! Neutrino
3014       part_charge( 5) =  1      ! Muon+
3015       part_charge( 6) = -1      ! Muon-
3016       part_charge( 7) =  0      ! Pion0
3017       part_charge( 8) =  1      ! Pion+
3018       part_charge( 9) = -1      ! Pion-
3019       part_charge(10) =  0      ! Kaon 0 long
3020       part_charge(11) =  1      ! Kaon+
3021       part_charge(12) = -1      ! Kaon-
3022       part_charge(13) =  0      ! Neutron
3023       part_charge(14) =  1      ! Proton
3024       part_charge(15) = -1      ! Antiproton
3025       part_charge(16) =  0      ! Kaon 0 short
3026       part_charge(17) =  0      ! Eta
3027       part_charge(18) =  0      ! Lambda
3028       part_charge(19) =  1      ! Sigma+
3029       part_charge(20) =  0      ! Sigma0
3030       part_charge(21) = -1      ! Sigma-
3031       part_charge(22) =  0      ! Xi 0
3032       part_charge(23) = -1      ! Xi -
3033       part_charge(24) = -1      ! Omega
3034       part_charge(25) =  0      ! Antineutron
3035       part_charge(26) =  0      ! Antilambda
3036       part_charge(27) = -1      ! Anti-Sigma -
3037       part_charge(28) =  0      ! Anti-Sigma 0
3038       part_charge(29) =  1      ! Anti-Sigma +
3039       part_charge(30) =  0      ! AntiXi 0
3040       part_charge(31) =  1      ! AntiXi +
3041       part_charge(32) =  1      ! Anti-Omega +
3042       part_charge(33) =  0
3043       part_charge(34) =  0
3044       part_charge(35) =  0
3045       part_charge(36) =  0
3046       part_charge(37) =  0
3047       part_charge(38) =  0
3048       part_charge(39) =  0
3049       part_charge(40) =  0
3050       part_charge(41) =  0
3051       part_charge(42) =  0
3052       part_charge(43) =  0
3053       part_charge(44) =  0
3054       part_charge(45) =  1      ! Deuteron
3055       part_charge(46) =  1      ! Triton
3056       part_charge(47) =  2      ! Alpha
3057       part_charge(48) =  0      ! Geantino (Fake particle)
3058       part_charge(49) =  2      ! He3
3059       part_charge(50) =  0      ! Cerenkov (Fake particle)
3060
3061 CCC   Set Particle Lifetimes:
3062
3063       part_lifetime( 1) = 1.0E+30       ! Gamma
3064       part_lifetime( 2) = 1.0E+30       ! Positron
3065       part_lifetime( 3) = 1.0E+30       ! Electron
3066       part_lifetime( 4) = 1.0E+30       ! Neutrino
3067       part_lifetime( 5) = 2.19703E-06   ! Muon+
3068       part_lifetime( 6) = 2.19703E-06   ! Muon-
3069       part_lifetime( 7) = 8.4E-17       ! Pion0
3070       part_lifetime( 8) = 2.603E-08     ! Pion+
3071       part_lifetime( 9) = 2.603E-08     ! Pion-
3072       part_lifetime(10) = 5.16E-08      ! Kaon 0 long
3073       part_lifetime(11) = 1.237E-08     ! Kaon+
3074       part_lifetime(12) = 1.237E-08     ! Kaon-
3075       part_lifetime(13) = 889.1         ! Neutron
3076       part_lifetime(14) = 1.0E+30       ! Proton
3077       part_lifetime(15) = 1.0E+30       ! Antiproton
3078       part_lifetime(16) = 8.922E-11     ! Kaon 0 short
3079       part_lifetime(17) = 5.53085E-19   ! Eta
3080       part_lifetime(18) = 2.632E-10     ! Lambda
3081       part_lifetime(19) = 7.99E-11      ! Sigma+
3082       part_lifetime(20) = 7.40E-20      ! Sigma0
3083       part_lifetime(21) = 1.479E-10     ! Sigma-
3084       part_lifetime(22) = 2.90E-10      ! Xi 0
3085       part_lifetime(23) = 1.639E-10     ! Xi -
3086       part_lifetime(24) = 8.22E-11      ! Omega
3087       part_lifetime(25) = 889.1         ! Antineutron
3088       part_lifetime(26) = 2.632E-10     ! Antilambda
3089       part_lifetime(27) = 7.99E-11      ! Anti-Sigma -
3090       part_lifetime(28) = 7.40E-20      ! Anti-Sigma 0
3091       part_lifetime(29) = 1.479E-10     ! Anti-Sigma +
3092       part_lifetime(30) = 2.90E-10      ! AntiXi 0
3093       part_lifetime(31) = 1.639E-10     ! AntiXi +
3094       part_lifetime(32) = 8.22E-11      ! Anti-Omega +
3095       part_lifetime(33) = 0.0
3096       part_lifetime(34) = 0.0
3097       part_lifetime(35) = 0.0
3098       part_lifetime(36) = 0.0
3099       part_lifetime(37) = 0.0
3100       part_lifetime(38) = 0.0
3101       part_lifetime(39) = 0.0
3102       part_lifetime(40) = 0.0
3103       part_lifetime(41) = 0.0
3104       part_lifetime(42) = 0.0
3105       part_lifetime(43) = 0.0
3106       part_lifetime(44) = 0.0
3107       part_lifetime(45) = 1.0E+30       ! Deuteron
3108       part_lifetime(46) = 1.0E+30       ! Triton
3109       part_lifetime(47) = 1.0E+30       ! Alpha
3110       part_lifetime(48) = 1.0E+30       ! Geantino (Fake particle)
3111       part_lifetime(49) = 1.0E+30       ! He3
3112       part_lifetime(50) = 1.0E+30       ! Cerenkov (Fake particle)
3113
3114       Return
3115       END
3116
3117 C----------------------------------------------------------------
3118
3119
3120       subroutine correl_model
3121       implicit none
3122
3123 CCC   This subroutine computes the requested 2-body model correlation
3124 CCC   function which is to be fitted by the track adjustment procedure.
3125 CCC   The model values are calculated on the requested fine and coarse
3126 CCC   mesh grid in momentum space.  The model values are computed at the
3127 CCC   mid point of each 1D bin or at the center of each 3D cell.  This
3128 CCC   could be refined at a later date to correspond to the integral of
3129 CCC   the model function over the bin width (cell volume) divided by the
3130 CCC   the bin width (cell volume).
3131
3132 C     The model includes the following options which are selected by the
3133 C     'switch*' parameters in common/parameters/:
3134 C
3135 C        switch_1d:          1D model as function of either Qinvar, Qtotal
3136 C                            or Q-vector
3137 C        switch_3d:          3D model as function of either the Bertsch-Pratt
3138 C                            side-out-long kinematics (but no cross term) or
3139 C                            the Yano-Koonin-Podgoretski perp-parallel-time
3140 C                            kinematics.
3141 C        switch_type:        Like and/or Unlike particles
3142 C        switch_coherence:   Purely incoherent source or a mixed incoherent-
3143 C                            coherent source.
3144 C        switch_coulomb:     Either (a) no Coulomb correction, (b) Gamow
3145 C                            factor, (c) NA35 parametrization, or (d) Pratt
3146 C                            Coulomb wave function integration for finite
3147 C                            size, spherical source.
3148 C        switch_fermi_bose:  Fermion or boson identical pairs.
3149
3150       Include 'common_parameters.inc'
3151       Include 'common_mesh.inc'
3152       Include 'common_correlations.inc'
3153
3154 CCC   Local Variable Type Declarations:
3155
3156       integer*4 i,j,k
3157
3158       real*4 R_1dsq, Rsidesq, Routsq, Rlongsq
3159       real*4 Rperpsq, Rparallelsq, R0sq
3160       real*4 sqrtlambda,fermi_bose_sign,coulomb_factor,coherence_fac
3161       real*4 q,q1,q2,q3
3162       real*4 b,b1,b2,b3
3163       real*4 massavg
3164
3165 CCC   Set Constants:
3166
3167       sqrtlambda  = sqrt(abs(lambda))
3168       R_1dsq      = ((R_1d     /hbc)**2)/2.0
3169       Rsidesq     = ((Rside    /hbc)**2)/2.0
3170       Routsq      = ((Rout     /hbc)**2)/2.0
3171       Rlongsq     = ((Rlong    /hbc)**2)/2.0
3172       Rperpsq     = ((Rperp    /hbc)**2)/2.0
3173       Rparallelsq = ((Rparallel/hbc)**2)/2.0
3174       R0sq        = ((R0       /hbc)**2)/2.0
3175
3176       fermi_bose_sign = float(switch_fermi_bose)
3177       coherence_fac   = switch_coherence*2.0*sqrtlambda*
3178      1                  (1.0 - sqrtlambda)
3179
3180 CCC   Determine average particle pair mass for Coulomb correction:
3181
3182       massavg = 0.14
3183       if(mass1.eq.0.0 .and. mass2.gt.0.0) massavg = mass2
3184       if(mass1.gt.0.0 .and. mass2.eq.0.0) massavg = mass1
3185       if(mass1.gt.0.0 .and. mass2.gt.0.0) massavg = 0.5*(mass1+mass2)
3186     
3187 CCC   Compute 1D correlation model arrays:
3188
3189       If(switch_1d .ge. 1) then
3190         If(n_1d_fine .gt. 0) then        ! Fill the 1D fine mesh bins
3191           q = -0.5*binsize_1d_fine
3192           do i = 1,n_1d_fine
3193             q = q + binsize_1d_fine
3194             b = exp(-q*q*R_1dsq)
3195
3196             if(switch_type.eq.1 .or. switch_type.eq.3) then
3197               c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3198      1                         *b*b + coherence_fac*b)
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_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3206             end if
3207
3208             if(switch_type.eq.2 .or. switch_type.eq.3) then
3209              c2mod_unlike_1d(i) = 1.0 
3210              if(switch_coulomb.eq.0) then
3211                coulomb_factor = 1.0
3212              else if(switch_coulomb.gt.0) then
3213                Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3214      1                      coulomb_factor)
3215              end if
3216              c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3217             end if
3218
3219           end do    ! End of 1D fine mesh filling do-loop
3220         end if      ! End of 1D fine mesh option
3221
3222         If(n_1d_coarse .gt. 0) then      ! Fill the 1D coarse mesh bins
3223           q = qmid_1d -0.5*binsize_1d_coarse
3224           do i = n_1d_fine + 1, n_1d_total
3225             q = q + binsize_1d_coarse
3226             b = exp(-q*q*R_1dsq)
3227
3228             if(switch_type.eq.1 .or. switch_type.eq.3) then
3229               c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3230      1                         *b*b + coherence_fac*b)
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_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3238             end if
3239
3240             if(switch_type.eq.2 .or. switch_type.eq.3) then
3241              c2mod_unlike_1d(i) = 1.0 
3242              if(switch_coulomb.eq.0) then
3243                coulomb_factor = 1.0
3244              else if(switch_coulomb.gt.0) then
3245                Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3246      1                      coulomb_factor)
3247              end if
3248              c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3249             end if
3250
3251           end do    ! End of 1D coarse mesh filling do-loop
3252         end if      ! End of 1D coarse mesh option
3253       end if        ! End of 1D option
3254
3255 CCC   Compute 3D correlation model arrays: