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
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.
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
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
33 C DESCRIPTION OF METHOD:
34 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:
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.
55 C (2) New one-body and two-body histograms, as well as the two-body
56 C correlation function are calculated.
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.
64 C (4) Steps 1-3 are repeated for each accepted particle in the
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.
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).
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/.
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.
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.
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.
131 C TWO-BODY REFERENCE HISTOGRAMS:
132 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
148 C ONE-BODY REFERENCE HISTOGRAMS:
149 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.
163 C TWO-BODY CORRELATION MODELS:
164 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:
172 C C2 = 1 + lambda*(b**2) + 2.0*sqrt(lambda)*[1 - sqrt(lambda)]*b
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:
179 C b = exp(- Q**2 * R**2 / 2)
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:
186 C b = exp[(- Qside**2 * Rside**2 - Qout**2 * Rout**2
187 C - Qlong**2 * Rlong**2)/2]
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):
193 C b = exp[(- Qperp**2 * Rperp**2 - Qparallel**2 * Rparallel**2
194 C - Qtime**2 * Rtime**2)/2]
197 C Qperp = transverse momentum difference
198 C Qparallel = Qlong = p_{1z} - p_{2z}
201 C The Coulomb correction may be omitted, or included in one of 3 ways:
203 C (1) Point source Gamow factor
204 C (2) Finite source NA35 model (see Eq.(5) in Z. Phys. C73, 443
207 C Coulomb correction = 1 + [G(q) - 1]*exp(-q/Q0)
209 C and G(q) is the Gamow factor and q is the 3-momentum
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.
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.
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.
231 C SUMMARY OF EXTERNAL FILES:
232 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 ----------------------------------------------------------------------------
251 C Source of Data for ALICE Application:
252 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 ----------------------------------------------------------------------------
271 C DESCRIPTION OF INPUT PARAMETERS AND SWITCHES (FILE: hbt_parameters.in):
272 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.
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
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.
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.
294 C switch_coherence = 0 for purely incoherent source (but can have
296 C = 1 for mixed incoherent and coherent source
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.
304 C switch_fermi_bose = 1 Boson pairs
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.
312 C print_full = 0 for standard, minimal output
313 C = 1 for full, comprehensive (large) output for
316 C print_sector_data = 0 std. sector occupancy data printed out
317 C = 1 to print sector occupancy and overflow info.
319 C Particle ID and Search Parameters:
320 C =================================
322 C n_pid_types = 1 or 2 only, # particle types to correlate
324 C pid(1), pid(2) = Geant-3 particle ID code numbers
326 C deltap = maximum range for random momentum shifts in
327 C GeV/c; px,py,pz independent; Def = 0.1 GeV/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.
335 C delchi = min % change in total chi-square for which
336 C the track adjustment iterations may stop,
339 C irand = initial random # seed, default = 12345
341 C Source Function Parameters:
342 C ==========================
344 C lambda = Chaoticity
346 C R_1d = Spherical source model radius (fm)
348 C Rside,Rout,Rlong = Non-spherical Bertsch-Pratt source model (fm)
350 C Rperp,Rparallel,R0= Non-spherical Yano-Koonin-Podgoretski source
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
360 C One-body pT, phi, eta Acceptance Bins:
361 C =====================================
363 C n_pt_bins, pt_min, pt_max = # pt bins, min/max pt accept. (GeV/c)
365 C n_phi_bins,phi_min,phi_max = # phi bins, min/max phi accept. (deg.)
367 C n_eta_bins,eta_min,eta_max = # eta bins, min/max eta accept.
369 C [NOTE: For each the maximum # of bins
372 C Two-body 1D and 3D Correlation Bins:
373 C ===================================
375 C n_1d_fine, binsize_1d_fine = # and size (GeV/c), 1D - fine mesh
377 C n_1d_coarse,binsize_1d_coarse = # and size (GeV/c), 1D - coarse mesh
379 C n_3d_fine, binsize_3d_fine = # and size (GeV/c), 3D - fine mesh
381 C n_3d_coarse,binsize_3d_coarse = # and size (GeV/c), 3D - coarse mesh
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
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.
394 C Momentum Space Track-Sector Cells:
395 C =================================
397 C n_px_bins,px_min,px_max = #, min,max px bins (GeV/c)
399 C n_py_bins,py_min,py_max = #, min,max py bins (GeV/c)
401 C n_pz_bins,pz_min,pz_max = #, min,max pz bins (GeV/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.]
409 C Relative Chi-Square Weights:
410 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
422 C FORMAT for hbt_singles_reference.hist:
423 C =====================================
425 C The output content for the one-body reference histograms is:
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
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)
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)
442 C FORMAT for hbt_pair_reference.hist:
443 C ==================================
445 C The output content for the two-body reference histograms is:
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
457 C The pair distributions (with one entry per line) are:
459 C 1D, like pairs: (href_like_1d(i),i=1,n_1d_total)
461 C 1D, unlike pairs: (href_unlike_1d(i),i=1,n_1d_total)
463 C 3D, like pairs, fine mesh: href_like_3d_fine(i,j,k) ; (i,(j,(k,...)))
465 C 3D, like pairs, coarse mesh: href_like_3d_coarse(i,j,k) ; (i,(j,(k,...)))
467 C 3D, unlike, fine mesh: href_unlike_3d_fine(i,j,k) ; (i,(j,(k,...)))
469 C 3D, unlike, coarse mesh: href_unlike_3d_coarse(i,j,k) ; (i,(j,(k,...)))
471 C*************************************************************************
472 C*************************************************************************
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'
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'
493 write(*,*) 'Input data in Fort'
497 write(*,*) ' PARAMETERS'
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
544 write(*,*) 'deltap',deltap
545 write(*,*) 'delchi',delchi
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
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
595 SUBROUTINE HBTPROCESSOR
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'
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'
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
615 CCC Initialize error code for ALICE application:
618 CCC Open Output Files:
620 open(unit=7,status='unknown',access='sequential',
621 1 name='hbt_log.out')
622 open(unit=8,status='unknown',access='sequential',
623 1 name='hbt_simulation.out')
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
635 CCC Read Input Controls and Parameters:
637 If(errorcode .eq. 1) Return
639 CCC Setup values and check input parameter ranges and consistency:
641 If(errorcode .eq. 1) Return
643 CCC Produce Basic Output File Header:
645 If(errorcode .eq. 1) Return
649 CCC Read Event Input file and fill flag files:
651 If(errorcode .eq. 1) Return
654 CCC Get the Reference Histograms and write out if new calculation:
656 If(errorcode .eq. 1) Return
658 If(errorcode .eq. 1) Return
662 CCC Compute the correlation model and print out:
665 If(errorcode .eq. 1) Return
668 CCC Carry out the Track Adjustment/Correlation Fitting Procedure:
672 CCC Final Output of Inclusive Quantities:
674 If(errorcode .eq. 1) Return
676 CCC Close Output Files:
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:')
692 C-------------------------------------------------------------------
695 subroutine initialize
698 CCC This subroutine sets all arrays and structures to zero:
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'
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'
712 CCC Local Variable Type Declarations:
724 trk_merge_flag(i) = 0
726 trk_start_vertex(i) = 0
727 trk_stop_vertex(i) = 0
728 trk_event_line(i) = 0
746 trk2_merge_flag(i) = 0
748 trk2_start_vertex(i) = 0
749 trk2_stop_vertex(i) = 0
750 trk2_event_line(i) = 0
765 stm_track_id(j,i) = 0
771 stm2_n_trk_sec(i) = 0
773 do j = 1,max_trk_sec2
774 stm2_track_id(j,i) = 0
782 part_lifetime(i) = 0.0
785 do i = 1,max_trk_save
792 hist_unlike_1d(i) = 0
794 htmp_unlike_1d(i) = 0
796 href_unlike_1d(i) = 0
798 hinc_unlike_1d(i) = 0
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
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
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
877 c2_coul_like(i) = 0.0
878 c2_coul_unlike(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
903 C---------------------------------------------------------------------
906 subroutine set_values
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.
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'
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'
927 CCC Local Variable Type Declarations:
929 integer*4 iphistep, ptmaxsteps, iptstep
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
936 CCC Check Input Controls:
938 if(ref_control .lt. 1 .or. ref_control .gt. 2) then
939 write(7,101) ref_control
944 if(switch_1d .lt. 0 .or. switch_1d .gt. 3) then
945 write(7,102) switch_1d
950 if(switch_3d .lt. 0 .or. switch_3d .gt. 2) then
951 write(7,103) switch_3d
956 if(switch_type .lt. 1 .or. switch_type .gt. 3) then
957 write(7,104) switch_type
962 if(switch_coherence .lt. 0 .or. switch_coherence .gt. 1) then
963 write(7,105) switch_coherence
968 if(switch_coulomb .lt. 0 .or. switch_coulomb .gt. 3) then
969 write(7,106) switch_coulomb
974 if(switch_fermi_bose.ne.-1 .and. switch_fermi_bose.ne.1) then
975 write(7,107) switch_fermi_bose
980 if(n_pid_types .lt. 1 .or. n_pid_types .gt. 2) then
981 write(7,108) n_pid_types
986 if(switch_type .ge. 2 .and. n_pid_types .eq. 1) then
987 write(7,109) switch_type, n_pid_types
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)
1000 if(pid(1).eq.0 .and. pid(2).eq.0) then
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
1010 write(7,1093) pid(1), pid(2)
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)
1022 if(trk_accep .le. 0.0) then
1023 write(7,10941) trk_accep
1029 CCC Check Input Parameters:
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
1036 CCC Check Coulomb source radius in range for Pratt type Coulomb correction.
1038 if(switch_coulomb .eq. 3 .and. (Q0 .lt. coulradmin .or.
1039 1 Q0 .gt. coulradmax)) then
1045 CCC Load the Pratt type Coulomb correction if this form is selected.
1047 if(switch_coulomb .eq. 3 .and. (Q0 .ge. coulradmin .and.
1048 1 Q0 .le. coulradmax)) then
1052 CCC Check and determine the one-body distribution's binning:
1054 if(n_pt_bins .lt. 1 .or. n_pt_bins .gt. max_h_1d) then
1055 write(7,110) n_pt_bins
1060 if(n_phi_bins .lt. 1 .or. n_phi_bins .gt. max_h_1d) then
1061 write(7,111) n_phi_bins
1066 if(n_eta_bins .lt. 1 .or. n_eta_bins .gt. max_h_1d) then
1067 write(7,112) n_eta_bins
1072 if(pt_min .gt. pt_max .or. pt_min .lt. 0.0) then
1073 write(7,113) pt_min, pt_max
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
1085 if(eta_min .gt. eta_max) then
1086 write(7,115) eta_min, eta_max
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)
1095 CCC Check and determine the two-body distribution's binning:
1097 n_1d_total = n_1d_fine + n_1d_coarse
1098 n_3d_total = n_3d_fine + n_3d_coarse - 1
1100 if(switch_1d .gt. 0) then
1101 if(n_1d_fine .lt. 1) then
1102 write(7,116) n_1d_fine
1107 if(n_1d_coarse .lt. 1) then
1108 write(7,117) n_1d_coarse
1113 if(n_1d_total .gt. max_h_1d) then
1114 write(7,118) n_1d_total
1119 qmid_1d = binsize_1d_fine *float(n_1d_fine)
1120 qmax_1d = binsize_1d_coarse*float(n_1d_coarse) + qmid_1d
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
1130 if(n_3d_coarse .lt. 1 .or. n_3d_coarse .gt. max_h_3d) then
1131 write(7,120) n_3d_coarse
1136 qmid_3d = binsize_3d_fine *float(n_3d_fine)
1137 qmax_3d = binsize_3d_coarse*float(n_3d_coarse)
1139 if(abs(qmid_3d - binsize_3d_coarse) .gt. 0.00001) then
1140 write(7,121) qmid_3d, binsize_3d_coarse
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
1150 if(n_3d_fine_project .lt. 1) then
1151 write(7,1212) n_3d_fine_project
1152 n_3d_fine_project = 1
1156 CCC Check and determine Track-Sector Parameters:
1158 if(n_px_bins .lt. 1) then
1159 write(7,122) n_px_bins
1164 if(n_py_bins .lt. 1) then
1165 write(7,123) n_py_bins
1170 if(n_pz_bins .lt. 1) then
1171 write(7,124) n_pz_bins
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
1183 if(n_sectors .gt. sec_maxlen2 .and. ref_control .eq. 2) then
1184 write(7,1251) n_sectors
1189 if(trk_maxlen .ne. trk2_maxlen .and. ref_control .eq. 2) then
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
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)
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.
1211 CCC Check the z-momentum components:
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)
1224 if(pz1 .lt. pz_min .or. pz2 .gt. pz_max) then
1225 write(7,126) pz1,pz_min,pz2,pz_max
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.
1236 CCC Divide the pt and phi acceptance ranges into 24 equal steps:
1238 pt_step = (pt_max - pt_min)/24.0
1239 phi_step = (phi_max - phi_min)/24.0
1244 phi1 = phi_min - phi_step
1246 phi1 = phi1 + phi_step
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
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
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
1265 if(pxstepmin .lt. px_min .or. pxstepmax .gt. px_max) then
1266 write(7,127) pxstepmin,px_min,pxstepmax,px_max
1271 if(pystepmin .lt. py_min .or. pystepmax .gt. py_max) then
1272 write(7,128) pystepmin,py_min,pystepmax,py_max
1277 CCC Load Geant Particle Data:
1278 Call Hbtp_particle_prop
1280 CCC Check Requested Particle ID Numbers:
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)
1288 CCC Initialize Masses to 0.0
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
1299 mass1 = part_mass(pid(1))
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
1307 mass2 = part_mass(pid(2))
1309 else if(n_pid_types .eq. 2) then
1310 if(pid(1) .lt. 1 .or. pid(1) .gt. part_maxlen) then
1315 mass1 = part_mass(pid(1))
1317 if(pid(2) .lt. 1 .or. pid(2) .gt. part_maxlen) then
1322 mass2 = part_mass(pid(2))
1326 CCC Set Math Constants:
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,
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')
1389 C----------------------------------------------------------------------
1392 subroutine read_data(mode)
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.
1400 C The type of input is controlled by the value of 'mode'
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.)
1406 C mode = 1, read the parameter and switches input file
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.
1412 C mode = 3, read the reference pair and one-body histograms
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.
1418 C mode = 5, same as mode=4, except the accepted tracks are
1419 C loaded into the 'trk2' data structure.
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
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'
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.
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
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 ---------------------------------------------------------------------------
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'
1470 Include 'common_track.inc'
1471 Include 'common_track2.inc'
1472 Include 'common_particle.inc'
1476 CCC Local Variable Type Declarations:
1478 real*4 px,py,pz,E,pt,phi,eta,mass
1479 real*4 acheck(10), function(20)
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
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
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)
1501 CHARACTER*100 CHROOT
1502 CHARACTER*100 FILNAM
1508 CCC Begin Data Input Options:
1510 C------------------------
1511 IF (mode.eq.1) THEN ! Read Input parameters from File#1
1512 C------------------------
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
1521 open(unit=1,type='old',access='sequential',
1522 1 name='hbt_parameters.in')
1524 CCC Read Control Switches: (See Main program listing for complete
1525 CCC description of input parameters)
1527 read(1,*) ref_control
1530 read(1,*) switch_type
1531 read(1,*) switch_coherence
1532 read(1,*) switch_coulomb
1533 read(1,*) switch_fermi_bose
1535 read(1,*) print_full
1536 read(1,*) print_sector_data
1538 CCC Read Parameters:
1540 read(1,*) n_pid_types
1541 read(1,*) pid(1),pid(2)
1547 CCC Read Source Parameters:
1551 read(1,*) Rside, Rout, Rlong
1552 read(1,*) Rperp, Rparallel, R0
1555 CCC Read one-body {pt,phi,eta} bins:
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
1561 CCC Read two-body 1D and 3D bins:
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
1569 CCC Read momentum space track sector bins in {px,py,pz}:
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
1575 CCC Relative Chi-Square weights for track adjustment fitting:
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
1587 End If ! ALICE Data I/O Option
1589 C-----------------------------
1590 ELSE IF (mode.eq.2) THEN
1591 C-----------------------------
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
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
1604 If(ALICE .eq. 0) Then
1605 open(unit=2,type='old',access='sequential',
1606 1 name='event_text.in')
1607 open(unit=3,status='unknown',access='sequential',
1608 1 name='event_line.flags')
1610 CCC Set Event Counter:
1613 30 read(2,10,err=20,end=25) evg_label
1615 if(evg_label .eq. event_line) then
1616 event_counter = event_counter + 1
1618 else if(evg_label .eq. vertex_line) then
1620 else if(evg_label .eq. track_line) then
1622 else if(evg_label .eq. gener_line) then
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')
1637 End If ! ALICE Data I/O Option
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
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.
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
1659 Call AliHbtp_GetNumberEvents(n_events)
1661 Call AliHbtp_SetActiveEventNumber(i)
1663 Call AliHbtp_GetNumberTracks(ntracks)
1665 Call AliHbtp_GetTrack(j,flag,px,py,pz,ge_pid)
1667 CCC Check if this track's particle ID is one to be used
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
1682 write(7,62) trk_maxlen, event_counter
1683 62 Format(5x,'#tracks exceeds trk_maxlen = ',
1684 1 I6,' for event#',I4)
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)
1694 C write(*,*) ' FFF: 2 calling PutTrack j = ',j
1695 Call AliHbtp_PutTrack(j,flag0,px,py,pz,ge_pid)
1700 Else If(ALICE .eq. 0) Then
1702 open(unit=2,type='old',access='sequential',
1703 1 name='event_text.in')
1704 open(unit=3,type='old',access='sequential',
1705 1 name='event_line.flags')
1706 open(unit=4,status='unknown',access='sequential',
1707 1 name='event_tracks.select')
1709 CCC Set Event Counter:
1712 40 read(3,11,err=45,end=50) flag
1714 event_counter = event_counter + 1
1721 else if(flag.eq.3) then
1722 read(2,41,err=46,end=50) ge_pid,px,py,pz,tid,start_v,
1724 41 Format(7x,I6,3(1x,G12.5),4(1x,I6))
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
1731 CCC Check if this track's particle ID is one to be used
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
1746 write(7,621) trk_maxlen, event_counter
1747 621 Format(5x,'#tracks exceeds trk_maxlen = ',
1748 1 I6,' for event#',I4)
1753 if(pidok.eq.1 .and. accepok.eq.1) then
1754 track_counter = track_counter + 1
1760 end if ! End Flag=3 options
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,
1767 46 write(7,61) event_counter
1768 61 Format(5x,'Read error in event_text.in (2nd pass) at event#',I5,
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.
1782 End If ! ALICE Data I/O Option
1784 C-----------------------------
1785 ELSE IF(mode.eq.3) THEN
1786 C-----------------------------
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.
1793 open(unit=9,type='old',access='sequential',
1794 1 name='hbt_pair_reference.hist')
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
1805 CCC Determine if the Input Reference pair histograms are compatible
1806 CCC with the present run parameters:
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
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
1834 if(ref_check .eq. 0) then
1836 70 Format(5x,'Reference Pair Histogram Parameters not consistent',
1837 1 ' with present run conditions - STOP')
1840 else if(ref_check .eq. 1) then
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)
1846 if(switch_type.eq.2 .or. switch_type.eq.3) then
1847 read(9,*) (href_unlike_1d(i),i=1,n_1d_total)
1849 end if ! End 1D input option
1851 if(switch_3d.gt.0) then
1852 if(switch_type.eq.1 .or. switch_type.eq.3) then
1854 if(n_3d_fine.gt.0) then
1858 read(9,*) href_like_3d_fine(i,j,k)
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)
1876 if(switch_type.eq.2 .or. switch_type.eq.3) then
1878 if(n_3d_fine.gt.0) then
1882 read(9,*) href_unlike_3d_fine(i,j,k)
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)
1900 end if ! End 3D input option
1901 end if ! End Reference Check OK/Not OK Option
1905 CCC Next read the one-body histograms for 1 or 2 particle ID types:
1907 open(unit=11,type='old',access='sequential',
1908 1 name='hbt_singles_reference.hist')
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
1916 CCC Determine if Reference one-body histograms are compatible with
1917 CCC the present run conditions.
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
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
1934 if(ref_check .eq. 0) then
1936 71 Format(5x,'Reference One-Body Histogram parameters not ',
1937 1 'consistent with current run - STOP')
1940 else if(ref_check .eq. 1) then
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)
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)
1954 end if ! End one-body reference histogram input
1958 C-----------------------------
1959 ELSE IF(mode.eq.4) THEN
1960 C-----------------------------
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'.
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.
1972 CCC Initialize accepted track counters for this new event:
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)
1978 If(ALICE .eq. 1) Then
1979 Call AliHbtp_GetNumberTracks(ntracks)
1981 Call AliHbtp_GetTrack(i,flag,px,py,pz,ge_pid)
1983 if(flag.eq.flag4) then
1984 track_counter = track_counter + 1
1986 if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
1987 track_counter_1 = track_counter_1 + 1
1990 if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
1991 track_counter_2 = track_counter_2 + 1
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
2017 n_part_1_trk = track_counter_1
2018 n_part_2_trk = track_counter_2
2019 n_part_tot_trk = track_counter
2021 Else If(ALICE .eq. 0) Then
2023 80 read(4,11,err=81,end=82) flag
2024 event_line_counter = event_line_counter + 1
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
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
2036 track_counter = track_counter + 1
2038 if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2039 track_counter_1 = track_counter_1 + 1
2042 if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2043 track_counter_2 = track_counter_2 + 1
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
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
2079 84 Format(5x,'Read error from file event_tracks.select for mode=4',
2083 85 Format(5x,'Read error from file event_text.in for mode=4',
2086 End If ! ALICE Data I/O Option
2088 C-----------------------------
2089 ELSE IF(mode.eq.5) THEN
2090 C-----------------------------
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'.
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.
2102 CCC Initialize accepted track counters for this new event:
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)
2108 If(ALICE .eq. 1) Then
2109 Call AliHbtp_GetNumberTracks(ntracks)
2111 Call AliHbtp_GetTrack(i,flag,px,py,pz,ge_pid)
2113 if(flag.eq.flag4) then
2114 track_counter = track_counter + 1
2116 if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2117 track_counter_1 = track_counter_1 + 1
2120 if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2121 track_counter_2 = track_counter_2 + 1
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
2147 n_part_1_trk2 = track_counter_1
2148 n_part_2_trk2 = track_counter_2
2149 n_part_tot_trk2 = track_counter
2151 Else If(ALICE.eq.0) Then
2153 90 read(4,11,err=91,end=92) flag
2154 event_line_counter = event_line_counter + 1
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
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
2166 track_counter = track_counter + 1
2168 if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2169 track_counter_1 = track_counter_1 + 1
2172 if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2173 track_counter_2 = track_counter_2 + 1
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
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
2209 94 Format(5x,'Read error from file event_tracks.select for mode=5',
2213 95 Format(5x,'Read error from file event_text.in for mode=5',
2217 End If ! ALICE Data I/O Option
2219 C-----------------------------
2220 ELSE IF(mode.eq.6) THEN
2221 C-----------------------------
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/.
2227 if(switch_coulomb.eq.3 .and. Q0.ge.coulradmin .and.
2228 1 Q0.le.coulradmax) then
2230 CCC Initially, read and interpolate like pair Coulomb corrections:
2233 If(ALICE .eq. 1) then
2235 CALL GETENVF('ALICE_ROOT',CHROOT)
2236 LNROOT = LNBLNK(CHROOT)
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*,'**********************************'
2251 FILNAM=CHROOT(1:LNROOT)//'/data/cpp_06.dat'
2252 open(unit=21,type='old',access='sequential',
2255 FILNAM=CHROOT(1:LNROOT)//'/data/cpp_08.dat'
2256 open(unit=22,type='old',access='sequential',
2259 FILNAM=CHROOT(1:LNROOT)//'/data/cpp_10.dat'
2260 open(unit=23,type='old',access='sequential',
2263 FILNAM=CHROOT(1:LNROOT)//'/data/cpp_12.dat'
2264 open(unit=24,type='old',access='sequential',
2267 FILNAM=CHROOT(1:LNROOT)//'/data/cpp_14.dat'
2268 open(unit=25,type='old',access='sequential',
2271 FILNAM=CHROOT(1:LNROOT)//'/data/cpp_16.dat'
2272 open(unit=26,type='old',access='sequential',
2275 FILNAM=CHROOT(1:LNROOT)//'/data/cpp_18.dat'
2276 open(unit=27,type='old',access='sequential',
2280 open(unit=21,type='old',access='sequential',
2281 1 name='cpp_06.dat')
2282 open(unit=22,type='old',access='sequential',
2283 1 name='cpp_08.dat')
2284 open(unit=23,type='old',access='sequential',
2285 1 name='cpp_10.dat')
2286 open(unit=24,type='old',access='sequential',
2287 1 name='cpp_12.dat')
2288 open(unit=25,type='old',access='sequential',
2289 1 name='cpp_14.dat')
2290 open(unit=26,type='old',access='sequential',
2291 1 name='cpp_16.dat')
2292 open(unit=27,type='old',access='sequential',
2293 1 name='cpp_18.dat')
2297 do i = 1,max_c2_coul
2298 do j = 1,ncoulradsteps
2299 read(20+j,*) q_coul(i), function(j)
2301 Call AliHbtp_interp(coulradmin,coulradmax,coulradstep,
2302 1 ncoulradsteps,function,20,Q0,c2_coul_like(i))
2313 CCC Next read and interpolate the unlike pair Coulomb corrections:
2315 If(ALICE .eq. 1) then
2316 FILNAM=CHROOT(1:LNROOT)//'/data/cpm_06.dat'
2317 open(unit=31,type='old',access='sequential',
2320 FILNAM=CHROOT(1:LNROOT)//'/data/cpm_08.dat'
2321 open(unit=32,type='old',access='sequential',
2324 FILNAM=CHROOT(1:LNROOT)//'/data/cpm_10.dat'
2325 open(unit=33,type='old',access='sequential',
2328 FILNAM=CHROOT(1:LNROOT)//'/data/cpm_12.dat'
2329 open(unit=34,type='old',access='sequential',
2332 FILNAM=CHROOT(1:LNROOT)//'/data/cpm_14.dat'
2333 open(unit=35,type='old',access='sequential',
2336 FILNAM=CHROOT(1:LNROOT)//'/data/cpm_16.dat'
2337 open(unit=36,type='old',access='sequential',
2340 FILNAM=CHROOT(1:LNROOT)//'/data/cpm_18.dat'
2341 open(unit=37,type='old',access='sequential',
2345 open(unit=31,type='old',access='sequential',
2346 1 name='cpm_06.dat')
2347 open(unit=32,type='old',access='sequential',
2348 1 name='cpm_08.dat')
2349 open(unit=33,type='old',access='sequential',
2350 1 name='cpm_10.dat')
2351 open(unit=34,type='old',access='sequential',
2352 1 name='cpm_12.dat')
2353 open(unit=35,type='old',access='sequential',
2354 1 name='cpm_14.dat')
2355 open(unit=36,type='old',access='sequential',
2356 1 name='cpm_16.dat')
2357 open(unit=37,type='old',access='sequential',
2358 1 name='cpm_18.dat')
2361 do i = 1,max_c2_coul
2362 do j = 1,ncoulradsteps
2363 read(30+j,*) q_coul(i), function(j)
2365 Call AliHbtp_interp(coulradmin,coulradmax,coulradstep,
2366 1 ncoulradsteps,function,20,Q0,c2_coul_unlike(i))
2377 CCC Convert the input q values which are in MeV/c, to GeV/c:
2379 do i = 1,max_c2_coul
2380 q_coul(i) = 0.001*q_coul(i)
2386 C----------------------------
2387 ELSE IF(mode.eq.7) THEN
2388 C----------------------------
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.
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.
2402 CCC Initialize accepted track counters for this new event:
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)
2408 If(ALICE .eq. 1) Then
2409 Call AliHbtp_GetNumberTracks(ntracks)
2411 Call AliHbtp_GetTrack(i,flag,px,py,pz,ge_pid)
2413 if(flag.eq.flag4) then
2414 track_counter = track_counter + 1
2416 if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2417 track_counter_1 = track_counter_1 + 1
2420 if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2421 track_counter_2 = track_counter_2 + 1
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
2447 n_part_1_trk = track_counter_1
2448 n_part_2_trk = track_counter_2
2449 n_part_tot_trk = track_counter
2451 Else If(ALICE .eq. 0) Then
2453 CCC Open temporary files:
2455 open(unit=12,status='unknown',access='sequential',
2456 1 name='event_text_aux.in')
2457 open(unit=14,status='unknown',access='sequential',
2458 1 name='event_tracks_aux.select')
2460 100 read(4,11,err=101,end=102) flag
2461 event_line_counter = event_line_counter + 1
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
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
2484 track_counter = track_counter + 1
2486 if(eg_pid.eq.pid(1) .and. pid(1).gt.0) then
2487 track_counter_1 = track_counter_1 + 1
2490 if(eg_pid.eq.pid(2) .and. pid(2).gt.0) then
2491 track_counter_2 = track_counter_2 + 1
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
2516 read(2,10,err=103,end=102) comment_event_label
2517 write(12,10) comment_event_label
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
2532 104 Format(5x,'Read error from file event_tracks.select for mode=7',
2536 105 Format(5x,'Read error from file event_text.in for mode=7',
2540 End If ! ALICE Data I/O Option
2542 C----------------------------
2543 ELSE IF(mode.eq.8) THEN
2544 C----------------------------
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.
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
2557 CCC Initialize accepted track counters:
2561 If(ALICE .eq. 1) Then
2562 Call AliHbtp_GetNumberTracks(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),
2579 Else If(ALICE .eq. 0) Then
2581 CCC Open temporary, auxiliary files:
2583 open(unit=12,type='old',access='sequential',
2584 1 name='event_text_aux.in')
2585 open(unit=14,type='old',access='sequential',
2586 1 name='event_tracks_aux.select')
2588 120 read(14,11,err=121,end=122) flag
2589 file10_line_counter = file10_line_counter + 1
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) ,
2617 7 trk_ge_pid(track_counter)
2618 841 Format('TRACK:',1x,I6,3(1x,G12.5),4(1x,I6))
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)
2631 read(12,10,err=123,end=122) comment_event_label
2632 write(10,10) comment_event_label
2635 if(flag .ne. 1) go to 120 ! Return to S.N. 120 and read next record
2637 122 Close(unit=12,status='delete')
2638 Close(unit=14,status='delete')
2641 124 Format(5x,'Read error from file event_tracks_aux.select',
2642 1 ' for mode = 8 - STOP')
2645 125 Format(5x,'Read error from file event_text_aux.in',
2646 1 ' for mode = 8 - STOP')
2649 End If ! ALICE Data I/O Option
2652 END IF ! End of read_data mode selection options
2658 C-----------------------------------------------------------------------
2661 subroutine getref_hist
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
2672 Include 'common_parameters.inc'
2673 Include 'common_mesh.inc'
2674 Include 'common_histograms.inc'
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'
2682 CCC Local Variable Type Declarations:
2684 integer*4 i,ipt,iphi,ieta,sign_toggle
2686 if(ref_control .eq. 1) then
2688 CCC read pair and one-body reference histograms:
2690 else if(ref_control .eq. 2) then
2692 CCC calculate the pair and one-body histograms:
2693 CCC Open event and flag files:
2695 If(ALICE .eq. 0) Then
2696 open(unit=2,type='old',access='sequential',
2697 1 name='event_text.in')
2698 open(unit=4,type='old',access='sequential',
2699 1 name='event_tracks.select')
2702 CCC Initialize counters:
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
2710 CCC Read event header lines (no tracks are in this part):
2711 If(ALICE .eq. 0) Then
2715 CCC Set toggle switch to alternate between loading event tracks into
2716 CCC table 'trk' and table 'trk2':
2719 CCC Start Event Loop:
2722 If(ALICE .eq. 1) Then
2723 Call AliHbtp_SetActiveEventNumber(i)
2725 if(sign_toggle .eq. 1) then ! Put tracks into 'trk'
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
2733 do ipt = 1,n_pt_bins
2734 href1_pt_1(ipt) = href1_pt_1(ipt) + hist1_pt_1(ipt)
2737 do iphi = 1,n_phi_bins
2738 href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2741 do ieta = 1,n_eta_bins
2742 href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2746 if(pid(2) .gt. 0) then
2747 Call histog1(1,0,2,pid(2),0.0,0.0,0.0)
2748 n_part_used_2_ref = n_part_used_2_ref + n_part_used_2_trk
2750 do ipt = 1,n_pt_bins
2751 href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2754 do iphi = 1,n_phi_bins
2755 href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2758 do ieta = 1,n_eta_bins
2759 href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2763 else if(sign_toggle .eq. (-1)) then ! Put tracks into 'trk2'
2766 Call stm_build(2,0,0)
2767 if(pid(1) .gt. 0) then
2768 Call histog1(4,0,1,pid(1),0.0,0.0,0.0)
2769 n_part_used_1_ref = n_part_used_1_ref +n_part_used_1_trk2
2771 do ipt = 1,n_pt_bins
2772 href1_pt_1(ipt) = href1_pt_1(ipt) + hist1_pt_1(ipt)
2775 do iphi = 1,n_phi_bins
2776 href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2779 do ieta = 1,n_eta_bins
2780 href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2784 if(pid(2) .gt. 0) then
2785 Call histog1(4,0,2,pid(2),0.0,0.0,0.0)
2786 n_part_used_2_ref = n_part_used_2_ref +n_part_used_2_trk2
2788 do ipt = 1,n_pt_bins
2789 href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2792 do iphi = 1,n_phi_bins
2793 href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2796 do ieta = 1,n_eta_bins
2797 href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2801 end if ! End read and load to trk or trk2 option
2803 sign_toggle = -sign_toggle
2805 if(i .gt. 1) then ! Compute 2-body reference histograms
2806 Call histog2(4,0,0,0,0,0.0,0.0,0.0,0.0)
2807 num_pairs_like_ref = num_pairs_like_ref
2808 1 + n_part_used_1_trk * n_part_used_1_trk2
2809 2 + n_part_used_2_trk * n_part_used_2_trk2
2810 num_pairs_unlike_ref = num_pairs_unlike_ref
2811 1 + n_part_used_1_trk * n_part_used_2_trk2
2812 2 + n_part_used_2_trk * n_part_used_1_trk2
2815 end do ! End of Event Loop
2817 CCC Write out the pair and one-body reference Histograms:
2818 Call write_data(2,0)
2820 If(ALICE .eq. 0) Then
2825 end if ! End Reference Histogram read/calculate option
2830 C----------------------------------------------------------------------
2833 subroutine AliHbtp_interp(rmin,rmax,rstep,nrsteps,function,
2837 CCC This routine interpolates the function values and puts the result
2838 CCC into 'answer'. It uses 2,3 or 4 mesh points which must be equally
2839 CCC spaced. The method uses the Lagrange interpolation formulas given
2840 CCC in Abramowitz and Stegun, ``Handbook of Mathematical Functions,''
2841 CCC (Dover Publications, New York, 1970); pages 878-879.
2843 CCC Definition of Variables in the Argument List:
2845 CCC rmin = lower limit of independent variable for input function
2846 CCC rmax = upper limit of independent variable for input function
2847 CCC rstep = step size of independent variable
2848 CCC nrsteps = (redundant) # of input steps
2849 CCC function(ndim) = Array of function values to be interpolated
2850 CCC ndim = array dimension size in calling program
2851 CCC r = coordinate value of independent variable where interpolation
2853 CCC answer = interpolated value
2855 CCC The algorithm will use the maximum number of points in the
2856 CCC interpolation, up to a maximum of 4
2858 CCC If the requested coordinate value, r, is out-of-range, then
2859 CCC 'answer' is returned with a 0.0 value.
2861 CCC Local Variable Type Declarations:
2863 integer*4 ndim, nrsteps, ik
2865 real*4 rmin,rmax,rstep,r,answer,rshift,p
2866 real*4 function(ndim),w1,w2,w3,w4
2870 if(abs(((rmax-rmin)/float(nrsteps-1))-rstep).gt.0.00001) then
2871 write(7,10) rmin,rmax,rstep,nrsteps
2872 10 Format(5x,'Interp mesh inconsistent:',3E12.5,I5,
2879 if(r .lt. rmin .or. r .gt. rmax) then
2880 write(7,11) rmin,rmax,r
2881 11 Format(5x,'Interp called with r out-of-range =',3E12.5)
2886 CCC Begin interpolation:
2888 if(nrsteps .eq. 2) then
2889 p = (r - rmin)/rstep
2890 answer = (1.0 - p)*function(1) + p*function(2)
2891 else if(nrsteps .eq. 3) then
2892 p = (r - (rmin + rstep))/rstep
2893 answer = 0.5*p*(p-1.0)*function(1) + (1.0 - p*p)
2894 1 *function(2) + 0.5*p*(p+1.0)*function(3)
2895 else if(nrsteps .ge. 4) then
2898 if(rshift .le. rstep) then
2900 p = (rshift - rstep)/rstep
2901 else if(rshift .ge. (rmax - rstep - rmin)) then
2903 p = (rshift - (rmax - rmin - 2.0*rstep))/rstep
2905 ik = int(rshift/rstep + 1.000001)
2906 if(ik .le. 1) ik = 2
2907 if(ik .ge. (nrsteps-1)) ik = nrsteps - 2
2908 p = (rshift - float(ik-1)*rstep)/rstep
2911 w1 = -p*(p-1.0)*(p-2.0)/6.0
2912 w2 = (p*p-1.0)*(p-2.0)/2.0
2913 w3 = -p*(p+1.0)*(p-2.0)/2.0
2914 w4 = p*(p*p-1.0)/6.0
2916 answer = w1*function(ik-1) + w2*function(ik)
2917 1 + w3*function(ik+1) + w4*function(ik+2)
2918 end if ! End # interplation points option
2923 C--------------------------------------------------------------------
2926 subroutine Hbtp_particle_prop
2929 CCC Fill particle properties table /particle/ with Geant 3 particle ID
2930 CCC numbers, charge (in units of |e|), mass in GeV/c**2 and lifetime
2931 CCC in seconds. See the Geant 3.15 Manual User's Guide, pages: CONS
2934 Include 'common_particle.inc'
2936 CCC Local Variable Type Declarations:
2940 do i = 1,part_maxlen
2944 CCC Set Particle Masses:
2946 part_mass( 1) = 0.0 ! Gamma
2947 part_mass( 2) = 0.00051099906 ! Positron
2948 part_mass( 3) = 0.00051099906 ! Electron
2949 part_mass( 4) = 0.0 ! Neutrino
2950 part_mass( 5) = 0.105658389 ! Muon+
2951 part_mass( 6) = 0.105658389 ! Muon-
2952 part_mass( 7) = 0.1349743 ! Pion0
2953 part_mass( 8) = 0.1395679 ! Pion+
2954 part_mass( 9) = 0.1395679 ! Pion-
2955 part_mass(10) = 0.497671 ! Kaon 0 long
2956 part_mass(11) = 0.493646 ! Kaon+
2957 part_mass(12) = 0.493646 ! Kaon-
2958 part_mass(13) = 0.93956563 ! Neutron
2959 part_mass(14) = 0.93827231 ! Proton
2960 part_mass(15) = 0.93827231 ! Antiproton
2961 part_mass(16) = 0.497671 ! Kaon 0 short
2962 part_mass(17) = 0.54745 ! Eta
2963 part_mass(18) = 1.11563 ! Lambda
2964 part_mass(19) = 1.18937 ! Sigma+
2965 part_mass(20) = 1.19255 ! Sigma0
2966 part_mass(21) = 1.197465 ! Sigma-
2967 part_mass(22) = 1.31485 ! Xi 0
2968 part_mass(23) = 1.32133 ! Xi -
2969 part_mass(24) = 1.67243 ! Omega
2970 part_mass(25) = 0.93956563 ! Antineutron
2971 part_mass(26) = 1.11563 ! Antilambda
2972 part_mass(27) = 1.18937 ! Anti-Sigma -
2973 part_mass(28) = 1.19255 ! Anti-Sigma 0
2974 part_mass(29) = 1.197465 ! Anti-Sigma +
2975 part_mass(30) = 1.31485 ! AntiXi 0
2976 part_mass(31) = 1.32133 ! AntiXi +
2977 part_mass(32) = 1.67243 ! Anti-Omega +
2990 part_mass(45) = 1.875613 ! Deuteron
2991 part_mass(46) = 2.80925 ! Triton
2992 part_mass(47) = 3.727417 ! Alpha
2993 part_mass(48) = 0.0 ! Geantino (Fake particle)
2994 part_mass(49) = 2.80923 ! He3
2995 part_mass(50) = 0.0 ! Cerenkov (Fake particle)
2997 CCC Set Particle Charge:
2999 part_charge( 1) = 0 ! Gamma
3000 part_charge( 2) = 1 ! Positron
3001 part_charge( 3) = -1 ! Electron
3002 part_charge( 4) = 0 ! Neutrino
3003 part_charge( 5) = 1 ! Muon+
3004 part_charge( 6) = -1 ! Muon-
3005 part_charge( 7) = 0 ! Pion0
3006 part_charge( 8) = 1 ! Pion+
3007 part_charge( 9) = -1 ! Pion-
3008 part_charge(10) = 0 ! Kaon 0 long
3009 part_charge(11) = 1 ! Kaon+
3010 part_charge(12) = -1 ! Kaon-
3011 part_charge(13) = 0 ! Neutron
3012 part_charge(14) = 1 ! Proton
3013 part_charge(15) = -1 ! Antiproton
3014 part_charge(16) = 0 ! Kaon 0 short
3015 part_charge(17) = 0 ! Eta
3016 part_charge(18) = 0 ! Lambda
3017 part_charge(19) = 1 ! Sigma+
3018 part_charge(20) = 0 ! Sigma0
3019 part_charge(21) = -1 ! Sigma-
3020 part_charge(22) = 0 ! Xi 0
3021 part_charge(23) = -1 ! Xi -
3022 part_charge(24) = -1 ! Omega
3023 part_charge(25) = 0 ! Antineutron
3024 part_charge(26) = 0 ! Antilambda
3025 part_charge(27) = -1 ! Anti-Sigma -
3026 part_charge(28) = 0 ! Anti-Sigma 0
3027 part_charge(29) = 1 ! Anti-Sigma +
3028 part_charge(30) = 0 ! AntiXi 0
3029 part_charge(31) = 1 ! AntiXi +
3030 part_charge(32) = 1 ! Anti-Omega +
3043 part_charge(45) = 1 ! Deuteron
3044 part_charge(46) = 1 ! Triton
3045 part_charge(47) = 2 ! Alpha
3046 part_charge(48) = 0 ! Geantino (Fake particle)
3047 part_charge(49) = 2 ! He3
3048 part_charge(50) = 0 ! Cerenkov (Fake particle)
3050 CCC Set Particle Lifetimes:
3052 part_lifetime( 1) = 1.0E+30 ! Gamma
3053 part_lifetime( 2) = 1.0E+30 ! Positron
3054 part_lifetime( 3) = 1.0E+30 ! Electron
3055 part_lifetime( 4) = 1.0E+30 ! Neutrino
3056 part_lifetime( 5) = 2.19703E-06 ! Muon+
3057 part_lifetime( 6) = 2.19703E-06 ! Muon-
3058 part_lifetime( 7) = 8.4E-17 ! Pion0
3059 part_lifetime( 8) = 2.603E-08 ! Pion+
3060 part_lifetime( 9) = 2.603E-08 ! Pion-
3061 part_lifetime(10) = 5.16E-08 ! Kaon 0 long
3062 part_lifetime(11) = 1.237E-08 ! Kaon+
3063 part_lifetime(12) = 1.237E-08 ! Kaon-
3064 part_lifetime(13) = 889.1 ! Neutron
3065 part_lifetime(14) = 1.0E+30 ! Proton
3066 part_lifetime(15) = 1.0E+30 ! Antiproton
3067 part_lifetime(16) = 8.922E-11 ! Kaon 0 short
3068 part_lifetime(17) = 5.53085E-19 ! Eta
3069 part_lifetime(18) = 2.632E-10 ! Lambda
3070 part_lifetime(19) = 7.99E-11 ! Sigma+
3071 part_lifetime(20) = 7.40E-20 ! Sigma0
3072 part_lifetime(21) = 1.479E-10 ! Sigma-
3073 part_lifetime(22) = 2.90E-10 ! Xi 0
3074 part_lifetime(23) = 1.639E-10 ! Xi -
3075 part_lifetime(24) = 8.22E-11 ! Omega
3076 part_lifetime(25) = 889.1 ! Antineutron
3077 part_lifetime(26) = 2.632E-10 ! Antilambda
3078 part_lifetime(27) = 7.99E-11 ! Anti-Sigma -
3079 part_lifetime(28) = 7.40E-20 ! Anti-Sigma 0
3080 part_lifetime(29) = 1.479E-10 ! Anti-Sigma +
3081 part_lifetime(30) = 2.90E-10 ! AntiXi 0
3082 part_lifetime(31) = 1.639E-10 ! AntiXi +
3083 part_lifetime(32) = 8.22E-11 ! Anti-Omega +
3084 part_lifetime(33) = 0.0
3085 part_lifetime(34) = 0.0
3086 part_lifetime(35) = 0.0
3087 part_lifetime(36) = 0.0
3088 part_lifetime(37) = 0.0
3089 part_lifetime(38) = 0.0
3090 part_lifetime(39) = 0.0
3091 part_lifetime(40) = 0.0
3092 part_lifetime(41) = 0.0
3093 part_lifetime(42) = 0.0
3094 part_lifetime(43) = 0.0
3095 part_lifetime(44) = 0.0
3096 part_lifetime(45) = 1.0E+30 ! Deuteron
3097 part_lifetime(46) = 1.0E+30 ! Triton
3098 part_lifetime(47) = 1.0E+30 ! Alpha
3099 part_lifetime(48) = 1.0E+30 ! Geantino (Fake particle)
3100 part_lifetime(49) = 1.0E+30 ! He3
3101 part_lifetime(50) = 1.0E+30 ! Cerenkov (Fake particle)
3106 C----------------------------------------------------------------
3109 subroutine correl_model
3112 CCC This subroutine computes the requested 2-body model correlation
3113 CCC function which is to be fitted by the track adjustment procedure.
3114 CCC The model values are calculated on the requested fine and coarse
3115 CCC mesh grid in momentum space. The model values are computed at the
3116 CCC mid point of each 1D bin or at the center of each 3D cell. This
3117 CCC could be refined at a later date to correspond to the integral of
3118 CCC the model function over the bin width (cell volume) divided by the
3119 CCC the bin width (cell volume).
3121 C The model includes the following options which are selected by the
3122 C 'switch*' parameters in common/parameters/:
3124 C switch_1d: 1D model as function of either Qinvar, Qtotal
3126 C switch_3d: 3D model as function of either the Bertsch-Pratt
3127 C side-out-long kinematics (but no cross term) or
3128 C the Yano-Koonin-Podgoretski perp-parallel-time
3130 C switch_type: Like and/or Unlike particles
3131 C switch_coherence: Purely incoherent source or a mixed incoherent-
3133 C switch_coulomb: Either (a) no Coulomb correction, (b) Gamow
3134 C factor, (c) NA35 parametrization, or (d) Pratt
3135 C Coulomb wave function integration for finite
3136 C size, spherical source.
3137 C switch_fermi_bose: Fermion or boson identical pairs.
3139 Include 'common_parameters.inc'
3140 Include 'common_mesh.inc'
3141 Include 'common_correlations.inc'
3143 CCC Local Variable Type Declarations:
3147 real*4 R_1dsq, Rsidesq, Routsq, Rlongsq
3148 real*4 Rperpsq, Rparallelsq, R0sq
3149 real*4 sqrtlambda,fermi_bose_sign,coulomb_factor,coherence_fac
3156 sqrtlambda = sqrt(abs(lambda))
3157 R_1dsq = ((R_1d /hbc)**2)/2.0
3158 Rsidesq = ((Rside /hbc)**2)/2.0
3159 Routsq = ((Rout /hbc)**2)/2.0
3160 Rlongsq = ((Rlong /hbc)**2)/2.0
3161 Rperpsq = ((Rperp /hbc)**2)/2.0
3162 Rparallelsq = ((Rparallel/hbc)**2)/2.0
3163 R0sq = ((R0 /hbc)**2)/2.0
3165 fermi_bose_sign = float(switch_fermi_bose)
3166 coherence_fac = switch_coherence*2.0*sqrtlambda*
3167 1 (1.0 - sqrtlambda)
3169 CCC Determine average particle pair mass for Coulomb correction:
3172 if(mass1.eq.0.0 .and. mass2.gt.0.0) massavg = mass2
3173 if(mass1.gt.0.0 .and. mass2.eq.0.0) massavg = mass1
3174 if(mass1.gt.0.0 .and. mass2.gt.0.0) massavg = 0.5*(mass1+mass2)
3176 CCC Compute 1D correlation model arrays:
3178 If(switch_1d .ge. 1) then
3179 If(n_1d_fine .gt. 0) then ! Fill the 1D fine mesh bins
3180 q = -0.5*binsize_1d_fine
3182 q = q + binsize_1d_fine
3183 b = exp(-q*q*R_1dsq)
3185 if(switch_type.eq.1 .or. switch_type.eq.3) then
3186 c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3187 1 *b*b + coherence_fac*b)
3188 if(switch_coulomb.eq.0) then
3189 coulomb_factor = 1.0
3190 else if(switch_coulomb.gt.0) then
3191 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3194 c2mod_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3197 if(switch_type.eq.2 .or. switch_type.eq.3) then
3198 c2mod_unlike_1d(i) = 1.0
3199 if(switch_coulomb.eq.0) then
3200 coulomb_factor = 1.0
3201 else if(switch_coulomb.gt.0) then
3202 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3205 c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3208 end do ! End of 1D fine mesh filling do-loop
3209 end if ! End of 1D fine mesh option
3211 If(n_1d_coarse .gt. 0) then ! Fill the 1D coarse mesh bins
3212 q = qmid_1d -0.5*binsize_1d_coarse
3213 do i = n_1d_fine + 1, n_1d_total
3214 q = q + binsize_1d_coarse
3215 b = exp(-q*q*R_1dsq)
3217 if(switch_type.eq.1 .or. switch_type.eq.3) then
3218 c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3219 1 *b*b + coherence_fac*b)
3220 if(switch_coulomb.eq.0) then
3221 coulomb_factor = 1.0
3222 else if(switch_coulomb.gt.0) then
3223 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3226 c2mod_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3229 if(switch_type.eq.2 .or. switch_type.eq.3) then
3230 c2mod_unlike_1d(i) = 1.0
3231 if(switch_coulomb.eq.0) then
3232 coulomb_factor = 1.0
3233 else if(switch_coulomb.gt.0) then
3234 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3237 c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3240 end do ! End of 1D coarse mesh filling do-loop
3241 end if ! End of 1D coarse mesh option
3242 end if ! End of 1D option
3244 CCC Compute 3D correlation model arrays:
3246 If(switch_3d .ge. 1) Then
3247 If(n_3d_fine .gt. 0) then ! Fill the 3D fine mesh bins
3248 q1 = -0.5*binsize_3d_fine
3250 q1 = q1 + binsize_3d_fine
3251 if(switch_3d.eq.1) b1=exp(-q1*q1*Rsidesq)
3252 if(switch_3d.eq.2) b1=exp(-q1*q1*Rperpsq)
3254 q2 = -0.5*binsize_3d_fine