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
3256 q2 = q2 + binsize_3d_fine
3257 if(switch_3d.eq.1) b2=exp(-q2*q2*Routsq)
3258 if(switch_3d.eq.2) b2=exp(-q2*q2*Rparallelsq)
3260 q3 = -0.5*binsize_3d_fine
3262 q3 = q3 + binsize_3d_fine
3263 if(switch_3d.eq.1) b3=exp(-q3*q3*Rlongsq)
3264 if(switch_3d.eq.2) b3=exp(-q3*q3*R0sq)
3267 if(switch_3d.eq.1) q = sqrt(q1*q1+q2*q2+q3*q3)
3268 if(switch_3d.eq.2) q = sqrt(q1*q1+q2*q2)
3270 if(switch_type.eq.1 .or. switch_type.eq.3) then
3271 c2mod_like_3d_fine(i,j,k) = 1.0 + fermi_bose_sign*(lambda
3272 1 *b*b + coherence_fac*b)
3273 if(switch_coulomb.eq.0) then
3274 coulomb_factor = 1.0
3275 else if(switch_coulomb.gt.0) then
3276 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3279 c2mod_like_3d_fine(i,j,k) =
3280 1 coulomb_factor*c2mod_like_3d_fine(i,j,k)
3283 if(switch_type.eq.2 .or. switch_type.eq.3) then
3284 c2mod_unlike_3d_fine(i,j,k) = 1.0
3285 if(switch_coulomb.eq.0) then
3286 coulomb_factor = 1.0
3287 else if(switch_coulomb.gt.0) then
3288 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3291 c2mod_unlike_3d_fine(i,j,k) =
3292 1 coulomb_factor*c2mod_unlike_3d_fine(i,j,k)
3297 end do ! End of 3D Fine Mesh Filling do-loops
3298 end if ! End 3D fine mesh option
3300 If(n_3d_coarse .gt. 0) then ! Fill the 3D coarse mesh bins
3301 q1 = -0.5*binsize_3d_coarse
3302 do i = 1,n_3d_coarse
3303 q1 = q1 + binsize_3d_coarse
3304 if(switch_3d.eq.1) b1=exp(-q1*q1*Rsidesq)
3305 if(switch_3d.eq.2) b1=exp(-q1*q1*Rperpsq)
3307 q2 = -0.5*binsize_3d_coarse
3308 do j = 1,n_3d_coarse
3309 q2 = q2 + binsize_3d_coarse
3310 if(switch_3d.eq.1) b2=exp(-q2*q2*Routsq)
3311 if(switch_3d.eq.2) b2=exp(-q2*q2*Rparallelsq)
3313 q3 = -0.5*binsize_3d_coarse
3314 do k = 1,n_3d_coarse
3315 q3 = q3 + binsize_3d_coarse
3316 if(switch_3d.eq.1) b3=exp(-q3*q3*Rlongsq)
3317 if(switch_3d.eq.2) b3=exp(-q3*q3*R0sq)
3320 if(switch_3d.eq.1) q = sqrt(q1*q1+q2*q2+q3*q3)
3321 if(switch_3d.eq.2) q = sqrt(q1*q1+q2*q2)
3323 if(switch_type.eq.1 .or. switch_type.eq.3) then
3324 c2mod_like_3d_coarse(i,j,k) = 1.0+fermi_bose_sign*(lambda
3325 1 *b*b + coherence_fac*b)
3326 if(switch_coulomb.eq.0) then
3327 coulomb_factor = 1.0
3328 else if(switch_coulomb.gt.0) then
3329 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3332 c2mod_like_3d_coarse(i,j,k) =
3333 1 coulomb_factor*c2mod_like_3d_coarse(i,j,k)
3336 if(switch_type.eq.2 .or. switch_type.eq.3) then
3337 c2mod_unlike_3d_coarse(i,j,k) = 1.0
3338 if(switch_coulomb.eq.0) then
3339 coulomb_factor = 1.0
3340 else if(switch_coulomb.gt.0) then
3341 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3344 c2mod_unlike_3d_coarse(i,j,k) =
3345 1 coulomb_factor*c2mod_unlike_3d_coarse(i,j,k)
3350 end do ! End of 3D Coarse Mesh Filling do-loops
3351 c2mod_like_3d_coarse(1,1,1) = 0.0
3352 c2mod_unlike_3d_coarse(1,1,1) = 0.0
3353 end if ! End of 3D Coarse Mesh Option
3355 End If ! End of 3D Option
3360 C----------------------------------------------------------------------
3363 subroutine coulomb(control,q,sign,mass,Q0,factor)
3366 CCC Compute Coulomb correction to the two-body correlation functions
3367 C for like and unlike charges for particles of the same mass.
3368 C Three methods are allowed:
3370 C If control = 1, Then use the gamow factor for point sources
3371 C If control = 2, Then use the NA35 finite source size empirical
3372 C correction factor from eq.(5) in Z. Phys. C73,
3374 C If control = 3, Then use the Pratt finite source size, numerically
3375 C integrated Coulomb correction factor with inter-
3378 C Other parameters in the argument list are:
3380 C q = 3-vector momentum difference for track pair, in GeV/c.
3381 C sign = algebraic sign of the charge product for the track pair.
3382 C mass = particle mass in GeV, it is assumed that both particles
3383 C have the same mass, e.g. pi+ and pi-, but not K+ and pi-.
3384 C Q0 = NA35 parameter in GeV/c if control = 2
3385 C = Source radius in fm if control = 3
3386 C factor = Multiplicative Coulomb correction result which is
3387 C calculated here and returned to the calling program.
3390 Include 'common_coulomb.inc'
3392 CCC Local Variable Type Declarations:
3394 integer*4 control, sign
3396 real*4 pi,q,mass,Q0,factor,alpha,eta,eta2pi
3398 parameter (pi = 3.141592654)
3399 parameter (alpha = 0.00729735)
3401 CCC Compute Gamow factor for control options 1 and 2:
3403 if(control.eq.1 .or. control.eq.2) then
3404 if(q .le. 0.001) then
3405 if(sign .gt. 0) gamow = 0.0
3406 if(sign .lt. 0) gamow = 86.0
3408 eta = sign*mass*alpha/q
3410 gamow = eta2pi/(exp(eta2pi) - 1.0)
3414 CCC Compute Coulomb Correction factor for options 1, 2 and 3:
3416 if(control .eq. 1) then
3418 else if(control .eq. 2) then
3419 factor = 1.0 + (gamow - 1.0)*exp(-q/Q0)
3420 else if(control .eq. 3) then
3422 if(q .le. q_coul(1)) then
3423 if(sign .gt. 0) factor = c2_coul_like(1)
3424 if(sign .lt. 0) factor = c2_coul_unlike(1)
3425 else if(q .ge. q_coul(max_c2_coul - 1)) then
3426 if(sign .gt. 0) factor = c2_coul_like(max_c2_coul - 1)
3427 if(sign .lt. 0) factor = c2_coul_unlike(max_c2_coul - 1)
3429 if(sign .gt. 0) then
3430 Call LAGRNG(q,q_coul,factor,c2_coul_like,
3431 1 max_c2_coul,1,5,max_c2_coul,1)
3432 else if(sign .lt. 0) then
3433 Call LAGRNG(q,q_coul,factor,c2_coul_unlike,
3434 1 max_c2_coul,1,5,max_c2_coul,1)
3437 end if ! END Coulomb correction evaluation, control selection opt.
3442 C---------------------------------------------------------------------
3445 SUBROUTINE LAGRNG (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
3446 IMPLICIT REAL*4(A-H,O-Z)
3448 C LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
3449 C ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
3450 C ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
3451 C VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
3452 C FUNCTIONS AT MAXARG VALUES.)
3453 C X =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
3454 C Y =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
3455 C NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
3456 C NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
3457 C NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
3459 DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
3461 C -----FIND X0, THE CLOSEST POINT TO X.
3465 10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
3466 IF ((NF-NI+1).EQ.2) GO TO 70
3468 IF (X.GT.ARG(NMID)) GO TO 20
3474 C ------ X IS ONE OF THE TABLULATED VALUES.
3476 30 IF (X.LE.ARG(NI)) GO TO 60
3485 C ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
3489 GO TO (110,100,90,80), NN
3491 IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
3495 IF ((N0+2).GT.NDIM) GO TO 110
3496 IF ((N0-2).LT.1) GO TO 100
3500 IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
3503 110 IF ((N0+1).LT.NDIM) GO TO 120
3505 C ------N0=NDIM, SPECIAL CASE.
3510 IF ((N0-1).LT.1) NUSED=2
3513 C ------AT LEAST 2 PTS LEFT.
3520 IF (NUSED.EQ.2) GO TO 140
3522 C ------AT LEAST 3 PTS.
3527 CM1=-Y0*Y1/Y0M1/YM11
3530 IF (NUSED.EQ.3) GO TO 160
3532 C ------AT LEAST 4 PTS
3541 C2=-YM1*Y0*Y1/YM12/Y02/Y12
3542 IF (NUSED.EQ.4) GO TO 180
3544 C ------AT LEAST 5 PTS.
3551 CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
3556 IF (NUSED.EQ.5) GO TO 200
3558 C ------AT LEAST 6 PTS.
3571 C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
3575 150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
3579 170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
3583 190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
3587 210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
3592 230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
3593 12*VAL(N,N0+2)+C3*VAL(N,N0+3)
3598 C-------------------------------------------------------------------
3601 subroutine Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,control)
3604 CCC Four-momentum kinematics conversion:
3606 C If control = 1, use input {px,py,pz,mass} to calculate
3608 C If control = 2, use input {pt,phi,eta,mass} to calculate
3611 C Units: Momentum are in GeV/c
3612 C Energy and mass are in GeV
3613 C Angles are in degrees
3615 CCC Local Variable Type Declarations:
3619 real*4 px,py,pz,E,pt,phi,eta,mass
3620 real*4 theta,pi,rad,pcut
3621 parameter (pi = 3.141592654)
3622 parameter (pcut = 0.000001)
3626 If(control .eq. 1) Then ! Use {px,py,pz,mass} --> {E,pt,phi,eta}
3627 pt = sqrt(px*px + py*py)
3628 E = sqrt(pt*pt + pz*pz + mass*mass)
3630 CCC Compute azimuthal angle phi; treat pt = 0.0 and py = 0.0 cases
3633 if(pt .le. pcut) then
3635 else if(pt.gt.pcut.and.abs(py).le.pcut.and.px.lt.0.0) then
3640 if(phi .lt. 0.0) phi = phi + 2.0*pi
3643 CCC Compute pseudorapidity:
3645 if(pt.le.pcut .and. abs(pz).le.pcut) then
3647 else if(pt.le.pcut .and. abs(pz).gt.pcut) then
3648 eta = 0.5*log((E+pz)/(E-pz)) ! Use beam rapidity
3650 theta = atan2(pt,pz)
3651 eta = -log(tan(theta/2.0))
3654 Else If(control .eq. 2) Then ! Use {pt,phi,eta,mass} --> {E,px,py,pz}
3656 px = pt*cos(phi/rad)
3657 py = pt*sin(phi/rad)
3658 if(abs(eta) .le. pcut) then
3661 theta = 2.0*atan(exp(-eta))
3665 E = sqrt(pt*pt + pz*pz + mass*mass)
3667 End If ! End control options
3672 C----------------------------------------------------------------------
3675 subroutine qdiff(px1,py1,pz1,E1,px2,py2,pz2,E2,qinvar2,qtotal2,
3676 1 qvector2,qside2,qout2,qlong2,qperp2,qtime2)
3679 CCC This subroutine computes the various relative momenta for given
3680 CCC input 4-momentum for particles 1 and 2. The subr: returns the
3681 CCC square of the momentum. All energy and momenta are in GeV.
3683 C Input 4-momentum for particle 1: {px1,py1,pz1,E1}
3684 C Input 4-momentum for particle 2: {px2,py2,pz2,E2}
3686 CCC Computed Momentum Differences are the following:
3688 C qinvar2 = Q-invariant**2
3689 C qtotal2 = Q-Total**2 (space**2 + time**2)
3690 C qvector2 = 3-momentum vector difference squared
3691 C qside2 = q-side**2 of Bertsch-Pratt 3D source models
3692 C qout2 = q-out**2 of Bertsch-Pratt 3D source models
3693 C qlong2 = q-long**2 of Bertsch-Pratt 3D source models
3694 C qperp2 = q-perpendicular**2 of YKP 3D source models
3695 C qparallel2= q-parallel**2 of YKP 3D source models
3696 C = qlong2 (not assigned a separate variable name)
3697 C qtime2 = q-time-like**2 of YKP 3D source models
3699 CCC Local Variable Type Declarations:
3701 real*4 px1,py1,pz1,E1,px2,py2,pz2,E2
3702 real*4 qinvar2,qtotal2,qvector2
3703 real*4 qside2,qout2,qlong2
3704 real*4 qperp2,qtime2
3705 real*4 px12sq,py12sq,pz12sq,E12sq
3707 px12sq = (px1 - px2)**2
3708 py12sq = (py1 - py2)**2
3709 pz12sq = (pz1 - pz2)**2
3710 E12sq = (E1 - E2)**2
3712 qvector2 = px12sq + py12sq + pz12sq
3713 qinvar2 = qvector2 - E12sq
3714 qtotal2 = qvector2 + E12sq
3717 qout2 = (px1*px1 - px2*px2 + py1*py1 - py2*py2)
3718 qout2 = qout2*qout2/((px1+px2)**2 + (py1+py2)**2)
3719 qperp2 = px12sq + py12sq
3720 qside2 = qperp2 - qout2
3722 if (qside2 .lt. 0) then
3723 C write(*,*) 'qside2 is less then 0', qside2
3724 C write(*,*) ' qperp2, qout2', qperp2, qout2
3725 C write(*,*) ' px1,py1,pz1,E1 ',px1,py1,pz1,E1
3726 C write(*,*) ' px2,py2,pz2,E2 ',px2,py2,pz2,E2
3732 C-----------------------------------------------------------------------
3735 subroutine mean_rms(a,ndim,npts,mean,rms)
3738 CCC Calculate the mean and standard deviation (rms) for input
3739 C distribution a() for npts number of values.
3740 C ndim = dimension of array a() in calling program.
3742 CCC Local Variable Type Declarations:
3744 integer*4 ndim, npts, i
3745 real*4 a(ndim), mean, rms, sum_mean, sum_rms
3747 if(npts .le. 0) then
3751 else if(npts .eq. 1) then
3759 sum_mean = sum_mean + a(i)
3761 mean = sum_mean/float(npts)
3764 sum_rms = sum_rms + (a(i) - mean)**2
3766 rms = sqrt((sum_rms)/float(npts - 1))
3772 C-----------------------------------------------------------------------
3775 subroutine tindex(mode,track_id)
3778 CCC This subroutine locates tracks in {px,py,pz} sectors
3779 C and sets the sector index numbers
3780 C in track table 'trk' and/or 'trk2', depending on the value of 'mode'.
3781 C If a track's momentum is out of the sector ranges, then the track
3782 C will be assigned to, and counted in the nearest sector cell on the
3785 C If mode = 1, apply to tracks in table 'trk'
3786 C If mode = 2, apply to tracks in table 'trk2'
3788 C If track_id = 0, then do this for all tracks
3789 C If track_id = i (where i.gt.0) then do this for track row i only
3791 Include 'common_parameters.inc'
3792 Include 'common_mesh.inc'
3794 Include 'common_track.inc'
3795 Include 'common_track2.inc'
3797 CCC Local Variable Type Declarations:
3799 integer*4 i,track_id,mode
3801 C-----------------------
3803 C-----------------------
3805 If(track_id .eq. 0) Then
3806 do i = 1,n_part_tot_trk
3807 trk_px_sec(i) = int(((trk_px(i) - px_min)/delpx)+1.00001)
3808 if(trk_px_sec(i) .lt.1) trk_px_sec(i) = 1
3809 if(trk_px_sec(i) .gt. n_px_bins) trk_px_sec(i) = n_px_bins
3810 trk_py_sec(i) = int(((trk_py(i) - py_min)/delpy)+1.00001)
3811 if(trk_py_sec(i) .lt.1) trk_py_sec(i) = 1
3812 if(trk_py_sec(i) .gt. n_py_bins) trk_py_sec(i) = n_py_bins
3813 trk_pz_sec(i) = int(((trk_pz(i) - pz_min)/delpz)+1.00001)
3814 if(trk_pz_sec(i) .lt.1) trk_pz_sec(i) = 1
3815 if(trk_pz_sec(i) .gt. n_pz_bins) trk_pz_sec(i) = n_pz_bins
3816 trk_sector(i) = trk_px_sec(i) + (trk_py_sec(i) - 1)*
3817 1 n_px_bins + (trk_pz_sec(i) - 1)*n_px_bins*n_py_bins
3820 Else If(track_id .gt. 0) Then
3822 trk_px_sec(i) = int(((trk_px(i) - px_min)/delpx)+1.00001)
3823 if(trk_px_sec(i) .lt.1) trk_px_sec(i) = 1
3824 if(trk_px_sec(i) .gt. n_px_bins) trk_px_sec(i) = n_px_bins
3825 trk_py_sec(i) = int(((trk_py(i) - py_min)/delpy)+1.00001)
3826 if(trk_py_sec(i) .lt.1) trk_py_sec(i) = 1
3827 if(trk_py_sec(i) .gt. n_py_bins) trk_py_sec(i) = n_py_bins
3828 trk_pz_sec(i) = int(((trk_pz(i) - pz_min)/delpz)+1.00001)
3829 if(trk_pz_sec(i) .lt.1) trk_pz_sec(i) = 1
3830 if(trk_pz_sec(i) .gt. n_pz_bins) trk_pz_sec(i) = n_pz_bins
3831 trk_sector(i) = trk_px_sec(i) + (trk_py_sec(i) - 1)*
3832 1 n_px_bins + (trk_pz_sec(i) - 1)*n_px_bins*n_py_bins
3835 C-----------------------------
3836 Else If (mode.eq.2) Then
3837 C-----------------------------
3839 If(track_id .eq. 0) Then
3840 do i = 1,n_part_tot_trk2
3841 trk2_px_sec(i) = int(((trk2_px(i) - px_min)/delpx)+1.00001)
3842 if(trk2_px_sec(i) .lt.1) trk2_px_sec(i) = 1
3843 if(trk2_px_sec(i) .gt. n_px_bins) trk2_px_sec(i) = n_px_bins
3844 trk2_py_sec(i) = int(((trk2_py(i) - py_min)/delpy)+1.00001)
3845 if(trk2_py_sec(i) .lt.1) trk2_py_sec(i) = 1
3846 if(trk2_py_sec(i) .gt. n_py_bins) trk2_py_sec(i) = n_py_bins
3847 trk2_pz_sec(i) = int(((trk2_pz(i) - pz_min)/delpz)+1.00001)
3848 if(trk2_pz_sec(i) .lt.1) trk2_pz_sec(i) = 1
3849 if(trk2_pz_sec(i) .gt. n_pz_bins) trk2_pz_sec(i) = n_pz_bins
3850 trk2_sector(i) = trk2_px_sec(i) + (trk2_py_sec(i) - 1)*
3851 1 n_px_bins + (trk2_pz_sec(i) - 1)*n_px_bins*n_py_bins
3854 Else If(track_id .gt. 0) Then
3856 trk2_px_sec(i) = int(((trk2_px(i) - px_min)/delpx)+1.00001)
3857 if(trk2_px_sec(i) .lt.1) trk2_px_sec(i) = 1
3858 if(trk2_px_sec(i) .gt. n_px_bins) trk2_px_sec(i) = n_px_bins
3859 trk2_py_sec(i) = int(((trk2_py(i) - py_min)/delpy)+1.00001)
3860 if(trk2_py_sec(i) .lt.1) trk2_py_sec(i) = 1
3861 if(trk2_py_sec(i) .gt. n_py_bins) trk2_py_sec(i) = n_py_bins
3862 trk2_pz_sec(i) = int(((trk2_pz(i) - pz_min)/delpz)+1.00001)
3863 if(trk2_pz_sec(i) .lt.1) trk2_pz_sec(i) = 1
3864 if(trk2_pz_sec(i) .gt. n_pz_bins) trk2_pz_sec(i) = n_pz_bins
3865 trk2_sector(i) = trk2_px_sec(i) + (trk2_py_sec(i) - 1)*
3866 1 n_px_bins + (trk2_pz_sec(i) - 1)*n_px_bins*n_py_bins
3876 C------------------------------------------------------------------------
3879 subroutine stm_build(mode,track_index,old_sector)
3882 CCC This subroutine fills or updates the track-sector information
3883 C table sec_trk_map or, for the reference calculations, it fills
3884 C sec_trk_map2. These track-sector tables contain the information
3885 C about the track occupancy, status, etc. for all of the {px,py,pz}
3890 C If track_index = 0, then fill information for all tracks in 'trk',
3892 C If track_index = i (where i.gt.0) then fill only the track-sector
3893 C information for track i in 'trk', into table 'stm'
3894 C Also for this case, if old_sector .ne. 0, then
3895 C remove the track # and ID information for this
3896 C old sector # from table stm
3898 C For Mode = 2: Fill information for all tracks in 'trk2' into table 'stm2'
3900 Include 'common_parameters.inc'
3901 Include 'common_mesh.inc'
3903 Include 'common_track.inc'
3904 Include 'common_track2.inc'
3905 Include 'common_sec_track.inc'
3906 Include 'common_sec_track2.inc'
3908 CCC Local Variable Type Declarations:
3910 integer*4 i,j,mode,track_index,old_sector,row
3911 integer*4 temp(max_trk_sec)
3913 C------------------------
3915 C------------------------
3917 If (track_index .eq. 0) Then
3920 stm_n_trk_sec(i) = 0
3921 do j = 1,max_trk_sec
3922 stm_track_id(j,i) = 0
3929 do i = 1,n_part_tot_trk
3930 if(trk_flag(i) .eq. 0) then
3932 stm_n_trk_sec(row) = stm_n_trk_sec(row) + 1
3933 if(stm_n_trk_sec(row) .le. max_trk_sec) then
3934 stm_track_id(stm_n_trk_sec(row),row) = trk_id(i)
3938 stm_n_trk_sec(row) = stm_n_trk_sec(row) - 1
3946 Else If (track_index .gt. 0) Then
3948 if(old_sector .ne. 0) then
3949 CCC Remove track from old sector:
3951 do i = 1,stm_n_trk_sec(old_sector)
3952 if(stm_track_id(i,old_sector) .ne. track_index) then
3954 temp(j) = stm_track_id(i,old_sector)
3957 stm_n_trk_sec(old_sector) = j
3958 do i = 1,max_trk_sec
3959 stm_track_id(i,old_sector) = 0
3961 do i = 1,stm_n_trk_sec(old_sector)
3962 stm_track_id(i,old_sector) = temp(i)
3965 CCC Update with new sector location of track:
3967 if(trk_flag(i) .eq. 0) then
3969 stm_n_trk_sec(row) = stm_n_trk_sec(row) + 1
3970 if(stm_n_trk_sec(row) .le. max_trk_sec) then
3971 stm_track_id(stm_n_trk_sec(row),row) = trk_id(i)
3975 stm_n_trk_sec(row) = stm_n_trk_sec(row) - 1
3983 C-----------------------------
3984 Else If (mode.eq.2) Then
3985 C-----------------------------
3987 If (track_index .eq. 0) Then
3988 do i = 1,sec_maxlen2
3990 stm2_n_trk_sec(i) = 0
3991 do j = 1,max_trk_sec2
3992 stm2_track_id(j,i) = 0
3999 do i = 1,n_part_tot_trk2
4000 if(trk2_flag(i) .eq. 0) then
4001 row = trk2_sector(i)
4002 stm2_n_trk_sec(row) = stm2_n_trk_sec(row) + 1
4003 if(stm2_n_trk_sec(row) .le. max_trk_sec2) then
4004 stm2_track_id(stm2_n_trk_sec(row),row) = trk2_id(i)
4008 stm2_n_trk_sec(row) = stm2_n_trk_sec(row) - 1
4018 End If ! End mode = 1,2 selection options
4024 C-----------------------------------------------------------------------
4027 subroutine sec_index(index,nbins,index_min,index_max)
4030 CCC Calculate track-sector neighboring bins and min->max range:
4032 CCC Local Variable Type Declarations:
4034 integer*4 index,nbins,index_min,index_max
4036 index_min = index - 1
4037 if(index_min .lt. 1) index_min = 1
4038 index_max = index + 1
4039 if(index_max .gt. nbins) index_max = nbins
4044 C-----------------------------------------------------------------------
4047 subroutine dist_range(mode,ntracks_out,ntracks_flagged)
4050 CCC Determine if tracks are out of acceptance range in pt, phi and eta,
4051 C and, if so, then set the 'out_flag' variable in the track table 'trk'
4053 CCC For Mode = 1, use track table 'trk'
4054 CCC For Mode = 2, use track table 'trk2'
4056 CCC Count the number of flagged tracks, i.e. trk(i).flag = 1, and "out
4057 C of acceptance range" tracks, i.e. trk(i).out_flag = 1, for both
4058 C particle ID types. Determine the number of tracks to use in the
4059 C correlation fit for each particle ID type.
4061 Include 'common_parameters.inc'
4062 Include 'common_mesh.inc'
4064 Include 'common_track.inc'
4065 Include 'common_track2.inc'
4067 CCC Local Variable Type Declarations:
4069 integer*4 i,mode,ntracks_out,ntracks_flagged
4071 C------------------------
4073 C------------------------
4081 do i = 1,n_part_tot_trk
4082 if(trk_flag(i) .eq. 0) then
4083 if(trk_pt(i) .lt. pt_min .or. trk_pt(i) .gt. pt_max)
4085 if(trk_phi(i).lt.phi_min .or. trk_phi(i).gt.phi_max)
4087 if(trk_eta(i).lt.eta_min .or. trk_eta(i).gt.eta_max)
4089 else if(trk_flag(i) .eq. 1) then
4090 ntracks_flagged = ntracks_flagged + 1
4095 do i = 1,n_part_tot_trk
4096 if(trk_out_flag(i) .eq. 1) ntracks_out = ntracks_out + 1
4099 n_part_used_1_trk = 0
4100 n_part_used_2_trk = 0
4101 do i = 1,n_part_tot_trk
4102 if(trk_flag(i) .eq. 0) then
4103 if(trk_ge_pid(i) .eq. pid(1)) then
4104 n_part_used_1_trk = n_part_used_1_trk + 1
4105 else if(trk_ge_pid(i) .eq. pid(2)) then
4106 n_part_used_2_trk = n_part_used_2_trk + 1
4111 C-----------------------------
4112 Else If (mode.eq.2) Then
4113 C-----------------------------
4115 do i = 1,trk2_maxlen
4116 trk2_out_flag(i) = 0
4121 do i = 1,n_part_tot_trk2
4122 if(trk2_flag(i) .eq. 0) then
4123 if(trk2_pt(i) .lt. pt_min .or. trk2_pt(i) .gt. pt_max)
4124 1 trk2_out_flag(i)=1
4125 if(trk2_phi(i).lt.phi_min .or. trk2_phi(i).gt.phi_max)
4126 1 trk2_out_flag(i)=1
4127 if(trk2_eta(i).lt.eta_min .or. trk2_eta(i).gt.eta_max)
4128 1 trk2_out_flag(i)=1
4129 else if(trk2_flag(i) .eq. 1) then
4130 ntracks_flagged = ntracks_flagged + 1
4135 do i = 1,n_part_tot_trk2
4136 if(trk2_out_flag(i) .eq. 1) ntracks_out = ntracks_out + 1
4139 n_part_used_1_trk2 = 0
4140 n_part_used_2_trk2 = 0
4141 do i = 1,n_part_tot_trk2
4142 if(trk2_flag(i) .eq. 0) then
4143 if(trk2_ge_pid(i) .eq. pid(1)) then
4144 n_part_used_1_trk2 = n_part_used_1_trk2 + 1
4145 else if(trk2_ge_pid(i) .eq. pid(2)) then
4146 n_part_used_2_trk2 = n_part_used_2_trk2 + 1
4152 End If ! End mode = 1,2 selection option
4158 C--------------------------------------------------------------------
4161 subroutine histog1(mode,itrack,pid_index,pidnum,pt_save,
4162 1 phi_save,eta_save)
4165 CCC This subroutine computes and/or updates the one-body histograms
4166 C according to the selected 'mode' and for the requested particle
4169 C Note: If the track momentum is out-of-range in {pt,phi,eta},
4170 C then it is ignored. The {pt,phi,eta} dependences for
4171 C the 1-dimensional histogramming are treated independently.
4172 C It is therefore possible for the sum of the number of
4173 C particles in the pt, phi and eta one-body, 1D histograms
4176 CCC Mode = 1, Fill the one-body histograms (hist1*) for selected
4177 C particle ID type, for the initial input distribution,
4178 C using the momenta in 'trk'
4180 C Mode = 2, Remove particle 'itrack' from temporary one-body hist-
4181 C ogram (htmp1*) for selected particle ID type, using
4182 C momentum values given by pt_save, phi_save, eta_save.
4184 C Mode = 3, Add particle 'itrack' to the temporary one-body hist-
4185 C ogram (htmp1*) for selected particle ID type, using
4186 C momentum values in track table 'trk'.
4188 C Mode = 4, Fill the one-body histograms (hist1*) for selected
4189 C particle ID type, for the initial input distribution,
4190 C using the momenta in 'trk2'
4192 C itrack = track index # for the removed or added track for mode =
4193 C 2 or 3, respectively.
4195 C pid_index = 1 or 2 for the first or second particle ID type, and
4196 C for filling/update either hist1*1 or hist1*2, and
4197 C similarly for htmp1*1 or htmp1*2
4199 C pidnum = Geant particle ID # for the track(s) to be filled or
4202 C pt_save = Removed track's pt value.
4204 C phi_save = Removed track's phi value.
4206 C eta_save = Removed track's eta value.
4208 Include 'common_parameters.inc'
4209 Include 'common_mesh.inc'
4210 Include 'common_histograms.inc'
4212 Include 'common_track.inc'
4213 Include 'common_track2.inc'
4215 CCC Local Variable Type Declarations:
4217 integer*4 mode,itrack,i,pid_index,pidnum,index
4218 integer*4 trk_counter,trk2_counter
4220 real*4 pt_save,phi_save,eta_save
4221 real*4 delpt,delphi,deleta
4223 C-------------------------
4225 C-------------------------
4227 CCC Fill one-body histograms for requested particle ID from table 'trk'
4228 CCC Initialize necessary arrays to zero:
4230 if(pid_index .eq. 1) then
4236 else if(pid_index .eq. 2) then
4246 do i = 1,n_part_tot_trk
4247 if(trk_ge_pid(i).eq.pidnum .and. trk_flag(i).eq.0) then
4248 trk_counter = trk_counter + 1
4249 delpt = trk_pt(i) - pt_min
4250 delphi = trk_phi(i) - phi_min
4251 deleta = trk_eta(i) - eta_min
4253 index = int((delpt/pt_bin_size) + 0.99999)
4254 if(index.ge.1 .and. index.le.max_h_1d) then
4255 if(pid_index.eq.1) hist1_pt_1(index) =
4256 1 hist1_pt_1(index) + 1
4257 if(pid_index.eq.2) hist1_pt_2(index) =
4258 1 hist1_pt_2(index) + 1
4261 index = int((delphi/phi_bin_size) + 0.99999)
4262 if(index.ge.1 .and. index.le.max_h_1d) then
4263 if(pid_index.eq.1) hist1_phi_1(index) =
4264 1 hist1_phi_1(index) + 1
4265 if(pid_index.eq.2) hist1_phi_2(index) =
4266 1 hist1_phi_2(index) + 1
4269 index = int((deleta/eta_bin_size) + 0.99999)
4270 if(index.ge.1 .and. index.le.max_h_1d) then
4271 if(pid_index.eq.1) hist1_eta_1(index) =
4272 1 hist1_eta_1(index) + 1
4273 if(pid_index.eq.2) hist1_eta_2(index) =
4274 1 hist1_eta_2(index) + 1
4280 if(pid_index .eq. 1) n_part_used_1_trk = trk_counter
4281 if(pid_index .eq. 2) n_part_used_2_trk = trk_counter
4283 C--------------------------------
4284 Else If (mode .eq. 2) Then
4285 C--------------------------------
4287 CCC Remove track # 'itrack' from fitting histograms in htmp1*,
4288 CCC use pt_save, phi_save, eta_save for the old momentum values
4289 CCC in order to determine which bins to decrement.
4291 if(trk_ge_pid(itrack).eq.pidnum.and.trk_flag(itrack).eq.0)then
4292 delpt = pt_save - pt_min
4293 delphi = phi_save - phi_min
4294 deleta = eta_save - eta_min
4296 index = int((delpt/pt_bin_size) + 0.99999)
4297 if(index.ge.1 .and. index.le.max_h_1d) then
4298 if(pid_index.eq.1) htmp1_pt_1(index) =
4299 1 htmp1_pt_1(index) - 1
4300 if(pid_index.eq.2) htmp1_pt_2(index) =
4301 1 htmp1_pt_2(index) - 1
4304 index = int((delphi/phi_bin_size) + 0.99999)
4305 if(index.ge.1 .and. index.le.max_h_1d) then
4306 if(pid_index.eq.1) htmp1_phi_1(index) =
4307 1 htmp1_phi_1(index) - 1
4308 if(pid_index.eq.2) htmp1_phi_2(index) =
4309 1 htmp1_phi_2(index) - 1
4312 index = int((deleta/eta_bin_size) + 0.99999)
4313 if(index.ge.1 .and. index.le.max_h_1d) then
4314 if(pid_index.eq.1) htmp1_eta_1(index) =
4315 1 htmp1_eta_1(index) - 1
4316 if(pid_index.eq.2) htmp1_eta_2(index) =
4317 1 htmp1_eta_2(index) - 1
4322 C--------------------------------
4323 Else If (mode .eq. 3) Then
4324 C--------------------------------
4326 CCC Add track # 'itrack' to fitting histograms in htmp1*,
4327 CCC use pt, phi and eta values in track table 'trk(itrack)'
4328 CCC for the new/added track position.
4330 if(trk_ge_pid(itrack).eq.pidnum.and.trk_flag(itrack).eq.0)then
4331 delpt = trk_pt(itrack) - pt_min
4332 delphi = trk_phi(itrack) - phi_min
4333 deleta = trk_eta(itrack) - eta_min
4335 index = int((delpt/pt_bin_size) + 0.99999)
4336 if(index.ge.1 .and. index.le.max_h_1d) then
4337 if(pid_index.eq.1) htmp1_pt_1(index) =
4338 1 htmp1_pt_1(index) + 1
4339 if(pid_index.eq.2) htmp1_pt_2(index) =
4340 1 htmp1_pt_2(index) + 1
4343 index = int((delphi/phi_bin_size) + 0.99999)
4344 if(index.ge.1 .and. index.le.max_h_1d) then
4345 if(pid_index.eq.1) htmp1_phi_1(index) =
4346 1 htmp1_phi_1(index) + 1
4347 if(pid_index.eq.2) htmp1_phi_2(index) =
4348 1 htmp1_phi_2(index) + 1
4351 index = int((deleta/eta_bin_size) + 0.99999)
4352 if(index.ge.1 .and. index.le.max_h_1d) then
4353 if(pid_index.eq.1) htmp1_eta_1(index) =
4354 1 htmp1_eta_1(index) + 1
4355 if(pid_index.eq.2) htmp1_eta_2(index) =
4356 1 htmp1_eta_2(index) + 1
4361 C------------------------------
4362 Else If (mode.eq.4) Then
4363 C------------------------------
4365 CCC Fill one-body histograms for requested particle ID from table 'trk2'
4366 CCC Initialize necessary arrays to zero:
4368 if(pid_index .eq. 1) then
4374 else if(pid_index .eq. 2) then
4384 do i = 1,n_part_tot_trk2
4385 if(trk2_ge_pid(i).eq.pidnum .and. trk2_flag(i).eq.0) then
4386 trk2_counter = trk2_counter + 1
4387 delpt = trk2_pt(i) - pt_min
4388 delphi = trk2_phi(i) - phi_min
4389 deleta = trk2_eta(i) - eta_min
4391 index = int((delpt/pt_bin_size) + 0.99999)
4392 if(index.ge.1 .and. index.le.max_h_1d) then
4393 if(pid_index.eq.1) hist1_pt_1(index) =
4394 1 hist1_pt_1(index) + 1
4395 if(pid_index.eq.2) hist1_pt_2(index) =
4396 1 hist1_pt_2(index) + 1
4399 index = int((delphi/phi_bin_size) + 0.99999)
4400 if(index.ge.1 .and. index.le.max_h_1d) then
4401 if(pid_index.eq.1) hist1_phi_1(index) =
4402 1 hist1_phi_1(index) + 1
4403 if(pid_index.eq.2) hist1_phi_2(index) =
4404 1 hist1_phi_2(index) + 1
4407 index = int((deleta/eta_bin_size) + 0.99999)
4408 if(index.ge.1 .and. index.le.max_h_1d) then
4409 if(pid_index.eq.1) hist1_eta_1(index) =
4410 1 hist1_eta_1(index) + 1
4411 if(pid_index.eq.2) hist1_eta_2(index) =
4412 1 hist1_eta_2(index) + 1
4418 if(pid_index .eq. 1) n_part_used_1_trk2 = trk2_counter
4419 if(pid_index .eq. 2) n_part_used_2_trk2 = trk2_counter
4422 End If ! End Mode = 1,2,3,4 Selection Options
4428 C-----------------------------------------------------------------------
4430 subroutine histog2(mode,itrack,px_sec_save,py_sec_save,
4431 1 pz_sec_save,px_save,py_save,pz_save,E_save)
4434 CCC This subroutine computes and/or updates the two-body histograms
4435 C according to the selected 'mode' and for the necessary particle
4438 C Mode = 1, Fill the two-body histograms (hist*) for all particles
4439 C in table 'trk' for like and unlike pairs, for 1D and/or
4440 C 3D fine and 3D coarse mesh bins.
4442 C Mode = 2, Remove all old track pairs for 'itrack' particle from
4443 C all htmp* histograms, for particles in table 'trk', for
4444 C like and unlike pairs, for 1D and/or 3D fine and 3D coarse
4445 C mesh bins; using the saved momentum and track values.
4447 C Mode = 3, Add all new track pairs for 'itrack' particle to
4448 C all htmp* histograms, for particles in table 'trk', for
4449 C like and unlike pairs, for 1D and/or 3D fine and 3D coarse
4450 C mesh bins; using the values in table 'trk(itrack).*'
4452 C Mode = 4, Fill and accumulate reference histograms (href*) for all
4453 C particle pairs from tables 'trk' and 'trk2', for like and
4454 C unlike pairs, for 1D and/or 3D fine and 3D coarse
4457 C itrack = single track index in table 'trk' for pairs to be removed
4458 C (mode = 2) or added (mode = 3).
4460 Include 'common_parameters.inc'
4461 Include 'common_mesh.inc'
4462 Include 'common_histograms.inc'
4464 Include 'common_track.inc'
4465 Include 'common_track2.inc'
4466 Include 'common_sec_track.inc'
4467 Include 'common_sec_track2.inc'
4469 CCC Local Variable Type Declarations:
4471 integer*4 mode,itrack,i,j,k,jx,jy,jz
4472 integer*4 jsec,trkj_sector,imin,imax,njtrks
4473 integer*4 index1,index2,index3
4474 integer*4 findex1,findex2,findex3
4475 integer*4 ipxmin,ipymin,ipzmin
4476 integer*4 ipxmax,ipymax,ipzmax
4477 integer*4 trki_pid,trkj_pid
4478 integer*4 px_sec_save,py_sec_save,pz_sec_save
4480 real*4 qinvar2,qtotal2,qvector2
4481 real*4 qside2, qout2, qlong2
4482 real*4 qperp2, qtime2
4483 real*4 qdiff1, qdiff2, qdiff3
4484 real*4 px_save,py_save,pz_save,E_save
4486 If (mode .eq. 1) Then ! Full hist* filling; initialize arrays to zero.
4490 hist_unlike_1d(i) = 0
4496 hist_like_3d_fine(i,j,k) = 0
4497 hist_unlike_3d_fine(i,j,k) = 0
4498 hist_like_3d_coarse(i,j,k) = 0
4499 hist_unlike_3d_coarse(i,j,k) = 0
4506 CCC Select # of particles to loop over for each 'mode':
4508 If (mode .eq. 1) Then
4510 imax = n_part_tot_trk
4511 Else If (mode .eq. 2 .or. mode .eq. 3) Then
4514 Else If (mode .eq. 4) Then
4516 imax = n_part_tot_trk
4519 C------------------------------------------------------
4520 CCC Begin Primary Loop over particles in Table 'trk':
4521 C------------------------------------------------------
4524 if(trk_flag(i) .eq. 0) then
4525 trki_pid = trk_ge_pid(i)
4527 Call sec_index(px_sec_save,n_px_bins,ipxmin,ipxmax)
4528 Call sec_index(py_sec_save,n_py_bins,ipymin,ipymax)
4529 Call sec_index(pz_sec_save,n_pz_bins,ipzmin,ipzmax)
4531 Call sec_index(trk_px_sec(i),n_px_bins,ipxmin,ipxmax)
4532 Call sec_index(trk_py_sec(i),n_py_bins,ipymin,ipymax)
4533 Call sec_index(trk_pz_sec(i),n_pz_bins,ipzmin,ipzmax)
4536 CCC Begin Loop over neighboring sectors for track # 'i':
4538 do jx = ipxmin,ipxmax
4539 do jy = ipymin,ipymax
4540 do jz = ipzmin,ipzmax
4541 trkj_sector = jx + (jy-1)*n_px_bins
4542 1 + (jz-1)*n_px_bins*n_py_bins
4544 if(mode.le.3) njtrks = stm_n_trk_sec(trkj_sector)
4545 if(mode.eq.4) njtrks = stm2_n_trk_sec(trkj_sector)
4546 if(njtrks .gt. 0) then
4548 CCC Begin Secondary Loop over particles in selected sectors in tables
4549 CCC 'trk' or 'trk2':
4551 if(mode.le.3) j = stm_track_id(jsec,trkj_sector)
4552 if(mode.eq.4) j = stm2_track_id(jsec,trkj_sector)
4553 if((mode.eq.1 .and. j.lt.i .and. trk_flag(j).eq.0)
4554 1 .or.(mode.eq.2 .and. j.ne.i .and. trk_flag(j).eq.0)
4555 2 .or.(mode.eq.3 .and. j.ne.i .and. trk_flag(j).eq.0)
4556 3 .or.(mode.eq.4 .and. trk2_flag(j).eq.0)) then
4558 CCC Obtain 1D and 3D relative momenta:
4560 if(mode.eq.1 .or. mode.eq.3) then
4561 trkj_pid = trk_ge_pid(j)
4562 Call qdiff(trk_px(i),trk_py(i),trk_pz(i),trk_E(i),
4563 1 trk_px(j),trk_py(j),trk_pz(j),trk_E(j),
4564 2 qinvar2,qtotal2,qvector2,qside2,qout2,qlong2,
4566 else if(mode.eq.2) then
4567 trkj_pid = trk_ge_pid(j)
4568 Call qdiff(px_save,py_save,pz_save,E_save,
4569 1 trk_px(j),trk_py(j),trk_pz(j),trk_E(j),
4570 2 qinvar2,qtotal2,qvector2,qside2,qout2,qlong2,
4572 else if(mode.eq.4) then
4573 trkj_pid = trk2_ge_pid(j)
4574 Call qdiff(trk_px(i),trk_py(i),trk_pz(i),trk_E(i),
4575 1 trk2_px(j),trk2_py(j),trk2_pz(j),trk2_E(j),
4576 2 qinvar2,qtotal2,qvector2,qside2,qout2,qlong2,
4580 C-----------------------------------------------
4581 CCC Fill and/or Update 1D two-body Histograms:
4582 C-----------------------------------------------
4584 if(switch_1d .gt. 0) then
4586 if(switch_1d .eq. 1) then
4587 qdiff1 = sqrt(qinvar2)
4588 else if(switch_1d .eq. 2) then
4589 qdiff1 = sqrt(qtotal2)
4590 else if(switch_1d .eq. 3) then
4591 qdiff1 = sqrt(qvector2)
4593 qdiff1 = sqrt(qvector2)
4596 if(qdiff1 .le. qmid_1d) then
4597 index1 = int((qdiff1/binsize_1d_fine)+0.99999)
4598 if(index1 .eq. 0) index1 = 1
4599 else if(qdiff1.gt.qmid_1d.and.qdiff1.le.qmax_1d) then
4600 index1 = int(((qdiff1-qmid_1d)/binsize_1d_coarse)
4602 if(index1.eq.0) index1 = 1
4603 index1 = index1 + n_1d_fine
4608 if(index1.ge.1.and.index1.le.n_1d_total) then
4609 if((trki_pid.eq.trkj_pid).and.(switch_type.eq.1
4610 1 .or. switch_type.eq.3)) then
4612 hist_like_1d(index1) = hist_like_1d(index1) + 1
4613 else if(mode.eq.2) then
4614 htmp_like_1d(index1) = htmp_like_1d(index1) - 1
4615 else if(mode.eq.3) then
4616 htmp_like_1d(index1) = htmp_like_1d(index1) + 1
4617 else if(mode.eq.4) then
4618 href_like_1d(index1) = href_like_1d(index1) + 1
4621 else if((trki_pid.ne.trkj_pid).and.(switch_type.eq.2
4622 1 .or. switch_type.eq.3)) then
4624 hist_unlike_1d(index1) = hist_unlike_1d(index1)+1
4625 else if(mode.eq.2) then
4626 htmp_unlike_1d(index1) = htmp_unlike_1d(index1)-1
4627 else if(mode.eq.3) then
4628 htmp_unlike_1d(index1) = htmp_unlike_1d(index1)+1
4629 else if(mode.eq.4) then
4630 href_unlike_1d(index1) = href_unlike_1d(index1)+1
4635 end if ! End 1D Histogram Fill and/or Update.
4637 C-----------------------------------------------
4638 CCC Fill and/or Update 3D Two-Body Histograms:
4639 C-----------------------------------------------
4641 if(switch_3d .gt. 0) then
4642 if(switch_3d .eq. 1) then
4643 qdiff1 = sqrt(qside2)
4644 qdiff2 = sqrt(qout2)
4645 qdiff3 = sqrt(qlong2)
4646 else if(switch_3d .eq. 2) then
4647 qdiff1 = sqrt(qperp2)
4648 qdiff2 = sqrt(qtime2)
4649 qdiff3 = sqrt(qlong2)
4651 qdiff1 = sqrt(qperp2)
4652 qdiff2 = sqrt(qtime2)
4653 qdiff3 = sqrt(qlong2)
4656 if(qdiff1 .le. qmid_3d) then
4657 findex1 = int((qdiff1/binsize_3d_fine)+0.99999)
4658 if(findex1 .eq. 0) findex1 = 1
4660 else if(qdiff1.gt.qmid_3d.and.qdiff1.le.qmax_3d) then
4661 index1 = int((qdiff1/binsize_3d_coarse)+0.99999)
4662 if(index1.eq.1) index1 = 2
4669 if(qdiff2 .le. qmid_3d) then
4670 findex2 = int((qdiff2/binsize_3d_fine)+0.99999)
4671 if(findex2 .eq. 0) findex2 = 1
4673 else if(qdiff2.gt.qmid_3d.and.qdiff2.le.qmax_3d) then
4674 index2 = int((qdiff2/binsize_3d_coarse)+0.99999)
4675 if(index2.eq.1) index2 = 2
4682 if(qdiff3 .le. qmid_3d) then
4683 findex3 = int((qdiff3/binsize_3d_fine)+0.99999)
4684 if(findex3 .eq. 0) findex3 = 1
4686 else if(qdiff3.gt.qmid_3d.and.qdiff3.le.qmax_3d) then
4687 index3 = int((qdiff3/binsize_3d_coarse)+0.99999)
4688 if(index3.eq.1) index3 = 2
4695 if((index1.ge.1.and.index1.le.n_3d_coarse).and.
4696 1 (index2.ge.1.and.index2.le.n_3d_coarse).and.
4697 2 (index3.ge.1.and.index3.le.n_3d_coarse)) then
4699 if((index1+index2+index3).eq.3) then
4701 if(findex1.ge.1.and.findex1.le.n_3d_fine.and.
4702 1 findex2.ge.1.and.findex2.le.n_3d_fine.and.
4703 2 findex3.ge.1.and.findex3.le.n_3d_fine) then
4705 if((trki_pid.eq.trkj_pid).and.(switch_type.eq.1
4706 1 .or. switch_type.eq.3)) then
4709 hist_like_3d_fine(findex1,findex2,findex3) =
4710 1 hist_like_3d_fine(findex1,findex2,findex3) +1
4711 else if(mode.eq.2) then
4712 htmp_like_3d_fine(findex1,findex2,findex3) =
4713 1 htmp_like_3d_fine(findex1,findex2,findex3) -1
4714 else if(mode.eq.3) then
4715 htmp_like_3d_fine(findex1,findex2,findex3) =
4716 1 htmp_like_3d_fine(findex1,findex2,findex3) +1
4717 else if(mode.eq.4) then
4718 href_like_3d_fine(findex1,findex2,findex3) =
4719 1 href_like_3d_fine(findex1,findex2,findex3) +1
4722 else if((trki_pid.ne.trkj_pid).and.(switch_type
4723 1 .eq.2 .or. switch_type.eq.3)) then
4726 hist_unlike_3d_fine(findex1,findex2,findex3) =
4727 1 hist_unlike_3d_fine(findex1,findex2,findex3) +1
4728 else if(mode.eq.2) then
4729 htmp_unlike_3d_fine(findex1,findex2,findex3) =
4730 1 htmp_unlike_3d_fine(findex1,findex2,findex3) -1
4731 else if(mode.eq.3) then
4732 htmp_unlike_3d_fine(findex1,findex2,findex3) =
4733 1 htmp_unlike_3d_fine(findex1,findex2,findex3) +1
4734 else if(mode.eq.4) then
4735 href_unlike_3d_fine(findex1,findex2,findex3) =
4736 1 href_unlike_3d_fine(findex1,findex2,findex3) +1
4742 else if((index1+index2+index3).gt.3) then
4744 if((trki_pid.eq.trkj_pid).and.(switch_type.eq.1
4745 1 .or. switch_type.eq.3)) then
4748 hist_like_3d_coarse(index1,index2,index3) =
4749 1 hist_like_3d_coarse(index1,index2,index3) +1
4750 else if(mode.eq.2) then
4751 htmp_like_3d_coarse(index1,index2,index3) =
4752 1 htmp_like_3d_coarse(index1,index2,index3) -1
4753 else if(mode.eq.3) then
4754 htmp_like_3d_coarse(index1,index2,index3) =
4755 1 htmp_like_3d_coarse(index1,index2,index3) +1
4756 else if(mode.eq.4) then
4757 href_like_3d_coarse(index1,index2,index3) =
4758 1 href_like_3d_coarse(index1,index2,index3) +1
4761 else if((trki_pid.ne.trkj_pid).and.(switch_type
4762 1 .eq.2 .or. switch_type.eq.3)) then
4765 hist_unlike_3d_coarse(index1,index2,index3) =
4766 1 hist_unlike_3d_coarse(index1,index2,index3) +1
4767 else if(mode.eq.2) then
4768 htmp_unlike_3d_coarse(index1,index2,index3) =
4769 1 htmp_unlike_3d_coarse(index1,index2,index3) -1
4770 else if(mode.eq.3) then
4771 htmp_unlike_3d_coarse(index1,index2,index3) =
4772 1 htmp_unlike_3d_coarse(index1,index2,index3) +1
4773 else if(mode.eq.4) then
4774 href_unlike_3d_coarse(index1,index2,index3) =
4775 1 href_unlike_3d_coarse(index1,index2,index3) +1
4779 end if ! End 3D Fine/Coarse Grid
4781 end if ! End 3D Histogram Filling and/or Update
4784 end do ! End Secondary Track Loop
4789 end do ! End Neighboring Sector Loop
4792 end do ! End Primary Track Loop
4797 C-----------------------------------------------------------------------
4800 subroutine correlation_fit
4803 CCC This subroutine carries out the track momentum adjustment
4804 CCC procedure in order to fit the requested model correlation
4805 CCC function and the input one-body {pt,phi,eta} distributions.
4807 CCC The method used is similar to the Metropolis method. Briefly:
4809 C 1. The accepted tracks for each event in the 'event_text.in'
4810 C input file are loaded into the 'trk' data structure table.
4811 C 2. The sector information, histograms, C2 and initial chi-square
4813 C 3. Each track momentum is randomly shifted within a specified
4814 C range, one track at a time, the histograms are updated, and
4815 C a new C2 and chi-square are computed.
4816 C 4. If the new track momentum is acceptable and if it results in a
4817 C smaller value of chi-square, then this shifted momentum is
4818 C retained, if not then the original momentum value is restored.
4819 C 5. This is done for all particles in the track table for the event.
4820 C 6. The entire process is repeated either until the maximum # of
4821 C iterations is reached, or until the chi-square improvement with
4822 C each iteration diminishes sufficiently.
4823 C 7. After completing the event loop, summary information is gathered
4824 C and inclusive correlation functions and one-body distributions
4827 Include 'common_parameters.inc'
4828 Include 'common_mesh.inc'
4829 Include 'common_histograms.inc'
4830 Include 'common_correlations.inc'
4831 Include 'common_event_summary.inc'
4833 Include 'common_track.inc'
4834 Include 'common_track2.inc'
4835 Include 'common_sec_track.inc'
4836 Include 'common_sec_track2.inc'
4837 Include 'common_particle.inc'
4839 CCC Local Variable Type Declarations:
4841 integer*4 i,j,k,ievent,niter,ntracks_out,nev
4842 integer*4 ntracks_flagged,track_status,pid_index
4843 integer*4 px_sec_save,py_sec_save,pz_sec_save
4845 real*4 px_save,py_save,pz_save,E_save
4846 real*4 pt_save,phi_save,eta_save,mass
4847 real*4 chisq_like_1d,chisq_unlike_1d
4848 real*4 chisq_like_3d_fine,chisq_unlike_3d_fine
4849 real*4 chisq_like_3d_coarse,chisq_unlike_3d_coarse
4850 real*4 chisq_hist1_1,chisq_hist1_2
4851 real*4 chisq_total,chisq_total_oldvalue,chisq_total_newvalue
4854 CCC Initialize counters:
4856 n_part_used_1_inc = 0
4857 n_part_used_2_inc = 0
4858 num_pairs_like_inc = 0
4859 num_pairs_unlike_inc = 0
4860 event_line_counter = 0
4861 file10_line_counter = 0
4863 CCC Open event input, track selection flags and event output files:
4865 If(ALICE .eq. 0) Then
4866 open(unit=2,type='old',access='sequential',
4867 1 name='event_text.in')
4868 open(unit=4,type='old',access='sequential',
4869 1 name='event_tracks.select')
4870 open(unit=10,status='unknown',access='sequential',
4871 1 name='event_hbt_text.out')
4872 CCC Read/Write event header from/to I/O event text files
4877 C-------------------------------------
4880 do ievent = 2, n_events + 1
4881 C-------------------------------------
4883 If(ALICE .eq. 1) Then
4884 Call AliHbtp_SetActiveEventNumber(ievent-1)
4885 C write(*,*) 'NEXT EVENT:', ievent
4888 if(n_part_tot_trk .gt. 0) then
4891 Call tindex(1,0) ! Fill initial track-sector info.
4892 Call stm_build(1,0,0) ! Fill initial sector info.
4893 Call dist_range(1,ntracks_out,ntracks_flagged)
4894 num_pairs_like = (n_part_used_1_trk*(n_part_used_1_trk-1))/2
4895 1 + (n_part_used_2_trk*(n_part_used_2_trk-1))/2
4896 num_pairs_unlike = n_part_used_1_trk*n_part_used_2_trk
4897 num_pairs_like_inc = num_pairs_like_inc + num_pairs_like
4898 num_pairs_unlike_inc = num_pairs_unlike_inc + num_pairs_unlike
4899 n_part_used_1_inc = n_part_used_1_inc + n_part_used_1_trk
4900 n_part_used_2_inc = n_part_used_2_inc + n_part_used_2_trk
4901 if(pid(1).gt.0) Call histog1(1,0,1,pid(1),0.,0.,0.)
4902 if(pid(2).gt.0) Call histog1(1,0,2,pid(2),0.,0.,0.)
4903 Call histog2(1,0,0,0,0,0.0,0.0,0.0,0.0)
4905 Call chisquare(1,chisq_like_1d,chisq_unlike_1d,
4906 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
4907 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
4908 3 chisq_hist1_1,chisq_hist1_2)
4909 chisq_total = chisq_wt_like_1d *chisq_like_1d
4910 1 + chisq_wt_unlike_1d *chisq_unlike_1d
4911 2 + chisq_wt_like_3d_fine *chisq_like_3d_fine
4912 3 + chisq_wt_unlike_3d_fine *chisq_unlike_3d_fine
4913 4 + chisq_wt_like_3d_coarse *chisq_like_3d_coarse
4914 5 + chisq_wt_unlike_3d_coarse *chisq_unlike_3d_coarse
4915 6 + chisq_wt_hist1_1 *chisq_hist1_1
4916 7 + chisq_wt_hist1_2 *chisq_hist1_2
4917 chisq_total_oldvalue = chisq_total
4922 1000 Continue ! Starting Point for Track Shift Iteration Loop:
4930 98 Format(5x,'** START NEXT EVENT **')
4931 99 Format(5x,'************************')
4932 write(8,100) ievent,niter,chisq_total
4933 100 Format(10x,'Event#+1,Iteration# and Chi-Sq = ',2I5,E11.4)
4935 IF(maxit .eq. 0) GO TO 1001 ! Option to compute correlations
4936 C ! for input events.
4938 C-------------------------------------
4939 C Begin Track Adjustment Loop:
4941 do i = 1,n_part_tot_trk
4942 C-------------------------------------
4944 if(trk_flag(i) .eq. 0) then
4946 CCC Save initial track parameters (those that might change):
4953 phi_save = trk_phi(i)
4954 eta_save = trk_eta(i)
4955 px_sec_save = trk_px_sec(i)
4956 py_sec_save = trk_py_sec(i)
4957 pz_sec_save = trk_pz_sec(i)
4958 old_sec_save = trk_sector(i)
4960 CCC Save the sector values that the track is initially located in:
4962 old_sec_ntrk = stm_n_trk_sec(trk_sector(i))
4963 old_sec_flag = stm_flag(trk_sector(i))
4964 do k = 1,stm_n_trk_sec(trk_sector(i))
4965 old_sec_trkid(k) = stm_track_id(k,trk_sector(i))
4968 CCC Determine the particle ID type:
4970 if(trk_ge_pid(i).eq.pid(1) .and. pid(1).gt.0) then
4972 else if(trk_ge_pid(i).eq.pid(2).and.pid(2).gt.0) then
4978 CCC Randomly shift track momentum vector and compute new kinematics:
4980 trk_px(i) = trk_px(i) + deltap*(2.0*hbtpran(irand) - 1.0)
4981 trk_py(i) = trk_py(i) + deltap*(2.0*hbtpran(irand) - 1.0)
4982 trk_pz(i) = trk_pz(i) + deltap*(2.0*hbtpran(irand) - 1.0)
4983 mass = part_mass(trk_ge_pid(i))
4984 Call Hbtp_kin(trk_px(i),trk_py(i),trk_pz(i),trk_E(i),
4985 1 trk_pt(i),trk_phi(i),trk_eta(i),mass,1)
4987 new_sec_save = trk_sector(i)
4989 CCC Determine if track has been shifted to a new sector, and if so,
4990 CCC whether this overfills this new sector. If all is well, then
4991 CCC update histograms. If not, then restore track parameters to their
4992 CCC initial values prior to shifting. Keep the status of track(i) in
4993 CCC 'track_status', where a value of 0 means the track is OK to use.
4995 CCC The Logical steps are the following:
4997 C IF(new track position is in same sector) THEN
4998 C o Remove old track position from htmp1*, htmp* using old saved values.
4999 C o Add new track position to htmp1*, htmp* using values in 'trk'
5000 C (Sector information is unchanged)
5001 C ELSE IF(new track position is in a different sector) THEN
5002 C IF(# tracks in new sector is still OK, with the new track) THEN
5003 C o Save values of new sector before trk was shifted into it.
5004 C o Remove old trk position from htmp1*, htmp* using old saved values
5005 C o Add new trk position to htmp1*, htmp* using values in trk
5006 C o Update sector information in stm
5007 C ELSE IF(# tracks in new sector becomes too many with new trk) THEN
5008 C o Restore track parameters to pre-shifted values
5009 C o Set track_status = 1, indicating the track could not be moved
5014 if(old_sec_save .eq. new_sec_save) then
5015 Call histog1(2,i,pid_index,pid(pid_index),pt_save,
5016 1 phi_save,eta_save)
5017 Call histog2(2,i,px_sec_save,py_sec_save,pz_sec_save,
5018 1 px_save,py_save,pz_save,E_save)
5020 Call histog1(3,i,pid_index,pid(pid_index),0.,0.,0.)
5021 Call histog2(3,i,0,0,0,0.0,0.0,0.0,0.0)
5023 else if(old_sec_save .ne. new_sec_save) then
5025 if(stm_n_trk_sec(new_sec_save) .lt. max_trk_sec) then
5026 new_sec_ntrk = stm_n_trk_sec(new_sec_save)
5027 new_sec_flag = stm_flag(new_sec_save)
5028 if(new_sec_ntrk .gt. 0) then
5029 do k = 1,new_sec_ntrk
5030 new_sec_trkid(k) = stm_track_id(k,new_sec_save)
5034 Call histog1(2,i,pid_index,pid(pid_index),pt_save,
5035 1 phi_save,eta_save)
5036 Call histog2(2,i,px_sec_save,py_sec_save,pz_sec_save,
5037 1 px_save,py_save,pz_save,E_save)
5039 Call histog1(3,i,pid_index,pid(pid_index),
5041 Call histog2(3,i,0,0,0,0.0,0.0,0.0,0.0)
5043 Call stm_build(1,i,old_sec_save)
5045 else if(stm_n_trk_sec(new_sec_save) .ge. max_trk_sec) then
5053 trk_phi(i) = phi_save
5054 trk_eta(i) = eta_save
5055 trk_px_sec(i) = px_sec_save
5056 trk_py_sec(i) = py_sec_save
5057 trk_pz_sec(i) = pz_sec_save
5058 trk_sector(i) = old_sec_save
5061 end if ! End Histogram and Sector Update
5063 CCC If the track was succesfully shifted then compute C2 and determine
5064 C if the chi-square decreases (improves) or increases. If it improves,
5065 C then store the new chi-square value and keep the 1- and 2-body
5066 C histograms in hist1* and hist*, repsectively. If chi-square
5067 C increases (worsens), then restore the track parameters to the
5068 C pre-shifted values, restore the histograms and if a new sector was
5069 C involved, then restore both the old and new sector values.
5071 C The Logical steps are the following:
5073 C IF(new track position is OK, (i.e. track_status = 0)) Then
5074 C o Compute C2 using htmp*
5075 C o Compute chi-square and save
5076 C IF(chi-square improves) Then
5077 C o Replace previous (best) chi-square with new value
5078 C o Update histograms, i.e. copy htmp1* -> hist1* and
5079 C copy htmp* -> hist*
5080 C ELSE IF(chi-square increases) Then
5081 C o Restore track parameters
5082 C o Restore histograms, i.e. copy hist1* -> htmp1* and
5083 C copy hist* -> htmp*
5084 C IF(new sector was used) Then
5085 C o Restore old sector values to pre-shifted values
5086 C o Restore new sector values to pre-shifted values
5091 If(track_status .eq.0) Then
5093 Call chisquare(2,chisq_like_1d,chisq_unlike_1d,
5094 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
5095 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
5096 3 chisq_hist1_1,chisq_hist1_2)
5097 chisq_total_newvalue =
5098 1 chisq_wt_like_1d *chisq_like_1d
5099 2 + chisq_wt_unlike_1d *chisq_unlike_1d
5100 3 + chisq_wt_like_3d_fine *chisq_like_3d_fine
5101 4 + chisq_wt_unlike_3d_fine *chisq_unlike_3d_fine
5102 5 + chisq_wt_like_3d_coarse *chisq_like_3d_coarse
5103 6 + chisq_wt_unlike_3d_coarse *chisq_unlike_3d_coarse
5104 7 + chisq_wt_hist1_1 *chisq_hist1_1
5105 8 + chisq_wt_hist1_2 *chisq_hist1_2
5107 if(chisq_total_newvalue .lt. chisq_total_oldvalue) then
5108 chisq_total_oldvalue = chisq_total_newvalue
5111 else if(chisq_total_newvalue.ge.chisq_total_oldvalue) then
5117 trk_phi(i) = phi_save
5118 trk_eta(i) = eta_save
5119 trk_px_sec(i) = px_sec_save
5120 trk_py_sec(i) = py_sec_save
5121 trk_pz_sec(i) = pz_sec_save
5122 trk_sector(i) = old_sec_save
5126 If(old_sec_save .ne. new_sec_save) then
5128 stm_n_trk_sec(old_sec_save) = old_sec_ntrk
5129 stm_flag(old_sec_save) = old_sec_flag
5130 do k = 1,max_trk_sec
5131 stm_track_id(k,old_sec_save) = 0
5133 do k = 1,old_sec_ntrk
5134 stm_track_id(k,old_sec_save) = old_sec_trkid(k)
5137 stm_n_trk_sec(new_sec_save) = new_sec_ntrk
5138 stm_flag(new_sec_save) = new_sec_flag
5139 do k = 1,max_trk_sec
5140 stm_track_id(k,new_sec_save) = 0
5142 do k = 1,new_sec_ntrk
5143 stm_track_id(k,new_sec_save) = new_sec_trkid(k)
5148 end if ! End Chi-Square Determination
5149 end if ! End valid tracks option
5150 end do ! End Track Shift Loop
5152 CCC Check chi-square for this iteration --
5153 C Best, current chi-square value is in 'chisq_total_oldvalue'
5154 C Chi-square value at the beginning of the iteration loop is in
5157 If(abs(200.0*(chisq_total_oldvalue - chisq_total)/
5158 1 (chisq_total_oldvalue + chisq_total)) .lt. delchi) then
5160 101 Format(/5x,'Chi-Sq reduced .lt. delchi % on last iteration',
5164 If (niter .gt. maxit) Then
5166 102 Format(/5x,'Max # Search Iterations Reached - Abort track ',
5170 chisq_total = chisq_total_oldvalue
5175 CCC Finished Track Adjustment Iteration Loop for event # 'ievent'
5177 if((ievent - 1) .le. max_events) then
5178 Call dist_range(1,ntracks_out,ntracks_flagged)
5179 num_iter(ievent-1) = float(niter)
5180 n_part_used_1_store(ievent-1) = float(n_part_used_1_trk)
5181 n_part_used_2_store(ievent-1) = float(n_part_used_2_trk)
5182 n_part_tot_store(ievent-1) = float(n_part_tot_trk)
5183 frac_trks_out(ievent-1)=float(ntracks_out)/
5184 1 float(n_part_tot_trk)
5185 frac_trks_flag(ievent-1) =
5186 1 float(ntracks_flagged)/float(n_part_tot_trk)
5190 Call chisquare(1,chisq_like_1d,chisq_unlike_1d,
5191 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
5192 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
5193 3 chisq_hist1_1,chisq_hist1_2)
5194 chisq_total = chisq_wt_like_1d *chisq_like_1d
5195 1 + chisq_wt_unlike_1d *chisq_unlike_1d
5196 2 + chisq_wt_like_3d_fine *chisq_like_3d_fine
5197 3 + chisq_wt_unlike_3d_fine *chisq_unlike_3d_fine
5198 4 + chisq_wt_like_3d_coarse *chisq_like_3d_coarse
5199 5 + chisq_wt_unlike_3d_coarse *chisq_unlike_3d_coarse
5200 6 + chisq_wt_hist1_1 *chisq_hist1_1
5201 7 + chisq_wt_hist1_2 *chisq_hist1_2
5203 if((ievent - 1) .le. max_events) then
5204 chisq_like_1d_store(ievent-1) = chisq_like_1d
5205 chisq_unlike_1d_store(ievent-1) = chisq_unlike_1d
5206 chisq_like_3d_fine_store(ievent-1) = chisq_like_3d_fine
5207 chisq_unlike_3d_fine_store(ievent-1) = chisq_unlike_3d_fine
5208 chisq_like_3d_coarse_store(ievent-1) = chisq_like_3d_coarse
5209 chisq_unlike_3d_coarse_store(ievent-1) =chisq_unlike_3d_coarse
5210 chisq_hist1_1_store(ievent-1) = chisq_hist1_1
5211 chisq_hist1_2_store(ievent-1) = chisq_hist1_2
5212 chisq_total_store(ievent-1) = chisq_total
5214 CCC Count # sectors with stm().flag = 1, indicating that too many
5215 C tracks were attempted to be loaded into that sector.
5217 num_sec_flagged_store(ievent-1) = 0.0
5219 if(stm_flag(k) .eq. 1) then
5220 num_sec_flagged_store(ievent-1) =
5221 1 num_sec_flagged_store(ievent-1) + 1.0
5228 if(print_full .eq. 1) Call write_data(5,ievent-1)
5230 end if ! End event-with-tracks processing.
5233 C-------------------------------
5234 end do ! End Event Loop
5235 C-------------------------------
5237 CCC Compute Correlation Functions for the Inclusive Histograms
5241 CCC Compute Mean and Std. dev of event monitor and summary quantities:
5243 if(n_events .le. max_events) then
5249 Call mean_rms(num_iter,nev,nev,niter_mean,niter_rms)
5250 Call mean_rms(n_part_used_1_store,nev,nev,npart1_mean,npart1_rms)
5251 Call mean_rms(n_part_used_2_store,nev,nev,npart2_mean,npart2_rms)
5252 Call mean_rms(n_part_tot_store,nev,nev,npart_tot_mean,
5254 Call mean_rms(num_sec_flagged_store,nev,nev,
5255 1 nsec_flag_mean,nsec_flag_rms)
5256 Call mean_rms(frac_trks_out,nev,nev,
5257 1 frac_trks_out_mean,frac_trks_out_rms)
5258 Call mean_rms(frac_trks_flag,nev,nev,
5259 1 frac_trks_flag_mean,frac_trks_flag_rms)
5260 Call mean_rms(chisq_like_1d_store,nev,nev,
5261 1 chi_l1d_mean,chi_l1d_rms)
5262 Call mean_rms(chisq_unlike_1d_store,nev,nev,
5263 1 chi_u1d_mean,chi_u1d_rms)
5264 Call mean_rms(chisq_like_3d_fine_store,nev,nev,
5265 1 chi_l3f_mean,chi_l3f_rms)
5266 Call mean_rms(chisq_unlike_3d_fine_store,nev,nev,
5267 1 chi_u3f_mean,chi_u3f_rms)
5268 Call mean_rms(chisq_like_3d_coarse_store,nev,nev,
5269 1 chi_l3c_mean,chi_l3c_rms)
5270 Call mean_rms(chisq_unlike_3d_coarse_store,nev,nev,
5271 1 chi_u3c_mean,chi_u3c_rms)
5272 Call mean_rms(chisq_hist1_1_store,nev,nev,
5273 1 chi_1_1_mean, chi_1_1_rms)
5274 Call mean_rms(chisq_hist1_2_store,nev,nev,
5275 1 chi_1_2_mean, chi_1_2_rms)
5276 Call mean_rms(chisq_total_store,nev,nev,
5277 1 chi_tot_mean, chi_tot_rms)
5279 If(ALICE .eq. 0) Then
5288 C------------------------------------------------------------------------
5291 subroutine hist1_copy(mode)
5294 CCC Copy 1-body histograms if:
5296 CCC mode = 1, then copy hist1* -> htmp1*
5297 CCC mode = 2, then copy htmp1* -> hist1*
5299 Include 'common_parameters.inc'
5300 Include 'common_mesh.inc'
5301 Include 'common_histograms.inc'
5303 CCC Local Variable Type Declarations:
5307 C---------------------------
5308 If(mode .eq. 1) Then ! Copy hist1* -> htmp1*
5309 C---------------------------
5311 if(pid(1) .gt. 0) then
5313 htmp1_pt_1(i) = hist1_pt_1(i)
5316 htmp1_phi_1(i) = hist1_phi_1(i)
5319 htmp1_eta_1(i) = hist1_eta_1(i)
5323 if(pid(2) .gt. 0) then
5325 htmp1_pt_2(i) = hist1_pt_2(i)
5328 htmp1_phi_2(i) = hist1_phi_2(i)
5331 htmp1_eta_2(i) = hist1_eta_2(i)
5335 C--------------------------------
5336 Else If (mode .eq. 2) Then ! Copy htmp1* -> hist1*
5337 C--------------------------------
5339 if(pid(1) .gt. 0) then
5341 hist1_pt_1(i) = htmp1_pt_1(i)
5344 hist1_phi_1(i) = htmp1_phi_1(i)
5347 hist1_eta_1(i) = htmp1_eta_1(i)
5351 if(pid(2) .gt. 0) then
5353 hist1_pt_2(i) = htmp1_pt_2(i)
5356 hist1_phi_2(i) = htmp1_phi_2(i)
5359 hist1_eta_2(i) = htmp1_eta_2(i)
5370 C----------------------------------------------------------------------
5372 subroutine hist2_copy(mode)
5375 CCC Copy 2-body histograms if:
5377 CCC mode = 1, then copy hist* -> htmp*
5378 CCC mode = 2, then copy htmp* -> hist*
5380 Include 'common_parameters.inc'
5381 Include 'common_mesh.inc'
5382 Include 'common_histograms.inc'
5384 CCC Local Variable Type Declarations:
5386 integer*4 mode, i,j,k
5388 C---------------------------
5389 If (mode .eq. 1) Then ! Copy hist* -> htmp*
5390 C---------------------------
5392 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
5393 if(switch_type.eq.1 .or. switch_type.eq.3) then
5395 htmp_like_1d(i) = hist_like_1d(i)
5398 if(switch_type.eq.2 .or. switch_type.eq.3) then
5400 htmp_unlike_1d(i) = hist_unlike_1d(i)
5403 end if ! End 1D histogram copy
5405 if(switch_3d.gt.0) then
5406 if(switch_type.eq.1 .or. switch_type.eq.3) then
5408 if(n_3d_fine .gt. 0) then
5412 htmp_like_3d_fine(i,j,k) = hist_like_3d_fine(i,j,k)
5418 if(n_3d_coarse .gt. 0) then
5419 do i = 1,n_3d_coarse
5420 do j = 1,n_3d_coarse
5421 do k = 1,n_3d_coarse
5422 htmp_like_3d_coarse(i,j,k) = hist_like_3d_coarse(i,j,k)
5430 if(switch_type.eq.2 .or. switch_type.eq.3) then
5432 if(n_3d_fine .gt. 0) then
5436 htmp_unlike_3d_fine(i,j,k) = hist_unlike_3d_fine(i,j,k)
5442 if(n_3d_coarse .gt. 0) then
5443 do i = 1,n_3d_coarse
5444 do j = 1,n_3d_coarse
5445 do k = 1,n_3d_coarse
5446 htmp_unlike_3d_coarse(i,j,k)=hist_unlike_3d_coarse(i,j,k)
5453 end if ! End 3D histogram copy
5455 C--------------------------------
5456 Else If (mode .eq. 2) Then ! Copy htmp* -> hist*
5457 C--------------------------------
5459 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
5460 if(switch_type.eq.1 .or. switch_type.eq.3) then
5462 hist_like_1d(i) = htmp_like_1d(i)
5465 if(switch_type.eq.2 .or. switch_type.eq.3) then
5467 hist_unlike_1d(i) = htmp_unlike_1d(i)
5470 end if ! End 1D histogram copy
5472 if(switch_3d.gt.0) then
5473 if(switch_type.eq.1 .or. switch_type.eq.3) then
5475 if(n_3d_fine .gt. 0) then
5479 hist_like_3d_fine(i,j,k) = htmp_like_3d_fine(i,j,k)
5485 if(n_3d_coarse .gt. 0) then
5486 do i = 1,n_3d_coarse
5487 do j = 1,n_3d_coarse
5488 do k = 1,n_3d_coarse
5489 hist_like_3d_coarse(i,j,k) = htmp_like_3d_coarse(i,j,k)
5497 if(switch_type.eq.2 .or. switch_type.eq.3) then
5499 if(n_3d_fine .gt. 0) then
5503 hist_unlike_3d_fine(i,j,k) = htmp_unlike_3d_fine(i,j,k)
5509 if(n_3d_coarse .gt. 0) then
5510 do i = 1,n_3d_coarse
5511 do j = 1,n_3d_coarse
5512 do k = 1,n_3d_coarse
5513 hist_unlike_3d_coarse(i,j,k)=htmp_unlike_3d_coarse(i,j,k)
5520 end if ! End 3D histogram copy
5523 End If ! End mode selection options
5529 C-----------------------------------------------------------------------
5532 subroutine hist1_incl_sum
5535 CCC Sum 1-body histograms for each event into inclusive totals, where
5536 CCC hinc1* = SUM[hist1*]
5538 Include 'common_parameters.inc'
5539 Include 'common_mesh.inc'
5540 Include 'common_histograms.inc'
5542 CCC Local Variable Type Declarations:
5546 if(pid(1) .gt. 0) then
5548 hinc1_pt_1(i) = hinc1_pt_1(i) + hist1_pt_1(i)
5551 hinc1_phi_1(i) = hinc1_phi_1(i) + hist1_phi_1(i)
5554 hinc1_eta_1(i) = hinc1_eta_1(i) + hist1_eta_1(i)
5558 if(pid(2) .gt. 0) then
5560 hinc1_pt_2(i) = hinc1_pt_2(i) + hist1_pt_2(i)
5563 hinc1_phi_2(i) = hinc1_phi_2(i) + hist1_phi_2(i)
5566 hinc1_eta_2(i) = hinc1_eta_2(i) + hist1_eta_2(i)
5574 C------------------------------------------------------------------------
5577 subroutine hist2_incl_sum
5580 CCC Sum 2-body histograms for each event into inclusive totals, where
5581 CCC hinc* = SUM[hist*]
5583 Include 'common_parameters.inc'
5584 Include 'common_mesh.inc'
5585 Include 'common_histograms.inc'
5587 CCC Local Variable Type Declarations:
5591 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
5592 if(switch_type.eq.1 .or. switch_type.eq.3) then
5594 hinc_like_1d(i) = hinc_like_1d(i) + hist_like_1d(i)
5597 if(switch_type.eq.2 .or. switch_type.eq.3) then
5599 hinc_unlike_1d(i) = hinc_unlike_1d(i) + hist_unlike_1d(i)
5602 end if ! End 1D Inclusive Histogram Sum
5604 if(switch_3d.gt.0) then
5605 if(switch_type.eq.1 .or. switch_type.eq.3) then
5607 if(n_3d_fine .gt. 0) then
5611 hinc_like_3d_fine(i,j,k) = hinc_like_3d_fine(i,j,k)
5612 1 + hist_like_3d_fine(i,j,k)
5618 if(n_3d_coarse .gt. 0) then
5619 do i = 1,n_3d_coarse
5620 do j = 1,n_3d_coarse
5621 do k = 1,n_3d_coarse
5622 hinc_like_3d_coarse(i,j,k) = hinc_like_3d_coarse(i,j,k)
5623 1 + hist_like_3d_coarse(i,j,k)
5631 if(switch_type.eq.2 .or. switch_type.eq.3) then
5633 if(n_3d_fine .gt. 0) then
5637 hinc_unlike_3d_fine(i,j,k) = hinc_unlike_3d_fine(i,j,k)
5638 1 + hist_unlike_3d_fine(i,j,k)
5644 if(n_3d_coarse .gt. 0) then
5645 do i = 1,n_3d_coarse
5646 do j = 1,n_3d_coarse
5647 do k = 1,n_3d_coarse
5648 hinc_unlike_3d_coarse(i,j,k) = hinc_unlike_3d_coarse(i,j,k)
5649 1 + hist_unlike_3d_coarse(i,j,k)
5656 end if ! End 3D Inclusive Histogram Sum
5661 C--------------------------------------------------------------------------
5664 subroutine correl_fit(mode)
5667 CCC This subroutine calculates the 2-body correlation function with
5668 CCC errors for the cases:
5670 CCC (1) 1D and/or 3D fine and coarse grid distributions
5671 CCC (2) like pairs and/or unlike pairs
5673 CCC Uses the signal and reference histograms. The input parameter
5674 CCC 'mode' selects which histograms to use.
5676 CCC Mode = 1, use hist*
5677 CCC Mode = 2, use htmp*
5678 CCC Mode = 3, use hinc*
5680 Include 'common_parameters.inc'
5681 Include 'common_mesh.inc'
5682 Include 'common_histograms.inc'
5683 Include 'common_correlations.inc'
5685 CCC Local Variable Type Declarations:
5687 integer*4 mode,i,j,k
5689 CCC Initialize correlation functions and error arrays to zero:
5692 c2fit_like_1d(i) = 0.0
5693 c2fit_unlike_1d(i) = 0.0
5694 c2err_like_1d(i) = 0.0
5695 c2err_unlike_1d(i) = 0.0
5701 c2fit_like_3d_fine(i,j,k) = 0.0
5702 c2fit_unlike_3d_fine(i,j,k) = 0.0
5703 c2fit_like_3d_coarse(i,j,k) = 0.0
5704 c2fit_unlike_3d_coarse(i,j,k) = 0.0
5705 c2err_like_3d_fine(i,j,k) = 0.0
5706 c2err_unlike_3d_fine(i,j,k) = 0.0
5707 c2err_like_3d_coarse(i,j,k) = 0.0
5708 c2err_unlike_3d_coarse(i,j,k) = 0.0
5713 CCC Compute 1D Correlation Functions and Errors:
5715 if(switch_1d .gt. 0) then
5716 if(switch_type.eq.1 .or. switch_type.eq.3) then
5718 if(mode .eq. 1) then
5719 Call c2_1d(hist_like_1d,href_like_1d,c2fit_like_1d,
5720 1 c2err_like_1d,max_h_1d,max_c2_1d,n_1d_total,
5721 2 num_pairs_like,num_pairs_like_ref)
5722 else if (mode .eq. 2) then
5723 Call c2_1d(htmp_like_1d,href_like_1d,c2fit_like_1d,
5724 1 c2err_like_1d,max_h_1d,max_c2_1d,n_1d_total,
5725 2 num_pairs_like,num_pairs_like_ref)
5726 else if (mode .eq. 3) then
5727 Call c2_1d(hinc_like_1d,href_like_1d,c2fit_like_1d,
5728 1 c2err_like_1d,max_h_1d,max_c2_1d,n_1d_total,
5729 2 num_pairs_like_inc,num_pairs_like_ref)
5734 if(switch_type.eq.2 .or. switch_type.eq.3) then
5736 if(mode .eq. 1) then
5737 Call c2_1d(hist_unlike_1d,href_unlike_1d,c2fit_unlike_1d,
5738 1 c2err_unlike_1d,max_h_1d,max_c2_1d,n_1d_total,
5739 2 num_pairs_unlike,num_pairs_unlike_ref)
5740 else if (mode .eq. 2) then
5741 Call c2_1d(htmp_unlike_1d,href_unlike_1d,c2fit_unlike_1d,
5742 1 c2err_unlike_1d,max_h_1d,max_c2_1d,n_1d_total,
5743 2 num_pairs_unlike,num_pairs_unlike_ref)
5744 else if (mode .eq. 3) then
5745 Call c2_1d(hinc_unlike_1d,href_unlike_1d,c2fit_unlike_1d,
5746 1 c2err_unlike_1d,max_h_1d,max_c2_1d,n_1d_total,
5747 2 num_pairs_unlike_inc,num_pairs_unlike_ref)
5750 end if ! End 1D correlations
5752 CCC Compute 3D Correlation Functions and Errors:
5754 if(switch_3d .gt. 0) then
5755 if(switch_type.eq.1 .or. switch_type.eq.3) then
5757 if(mode .eq. 1) then
5758 Call c2_3d(hist_like_3d_fine,href_like_3d_fine,
5759 1 c2fit_like_3d_fine,c2err_like_3d_fine,
5760 2 max_h_3d,max_c2_3d,n_3d_fine,
5761 3 num_pairs_like,num_pairs_like_ref)
5762 Call c2_3d(hist_like_3d_coarse,href_like_3d_coarse,
5763 1 c2fit_like_3d_coarse,c2err_like_3d_coarse,
5764 2 max_h_3d,max_c2_3d,n_3d_coarse,
5765 3 num_pairs_like,num_pairs_like_ref)
5766 else if(mode .eq. 2) then
5767 Call c2_3d(htmp_like_3d_fine,href_like_3d_fine,
5768 1 c2fit_like_3d_fine,c2err_like_3d_fine,
5769 2 max_h_3d,max_c2_3d,n_3d_fine,
5770 3 num_pairs_like,num_pairs_like_ref)
5771 Call c2_3d(htmp_like_3d_coarse,href_like_3d_coarse,
5772 1 c2fit_like_3d_coarse,c2err_like_3d_coarse,
5773 2 max_h_3d,max_c2_3d,n_3d_coarse,
5774 3 num_pairs_like,num_pairs_like_ref)
5775 else if(mode .eq. 3) then
5776 Call c2_3d(hinc_like_3d_fine,href_like_3d_fine,
5777 1 c2fit_like_3d_fine,c2err_like_3d_fine,
5778 2 max_h_3d,max_c2_3d,n_3d_fine,
5779 3 num_pairs_like_inc,num_pairs_like_ref)
5780 Call c2_3d(hinc_like_3d_coarse,href_like_3d_coarse,
5781 1 c2fit_like_3d_coarse,c2err_like_3d_coarse,
5782 2 max_h_3d,max_c2_3d,n_3d_coarse,
5783 3 num_pairs_like_inc,num_pairs_like_ref)
5788 if(switch_type.eq.2 .or. switch_type.eq.3) then
5790 if(mode .eq. 1) then
5791 Call c2_3d(hist_unlike_3d_fine,href_unlike_3d_fine,
5792 1 c2fit_unlike_3d_fine,c2err_unlike_3d_fine,
5793 2 max_h_3d,max_c2_3d,n_3d_fine,
5794 3 num_pairs_unlike,num_pairs_unlike_ref)
5795 Call c2_3d(hist_unlike_3d_coarse,href_unlike_3d_coarse,
5796 1 c2fit_unlike_3d_coarse,c2err_unlike_3d_coarse,
5797 2 max_h_3d,max_c2_3d,n_3d_coarse,
5798 3 num_pairs_unlike,num_pairs_unlike_ref)
5799 else if(mode .eq. 2) then
5800 Call c2_3d(htmp_unlike_3d_fine,href_unlike_3d_fine,
5801 1 c2fit_unlike_3d_fine,c2err_unlike_3d_fine,
5802 2 max_h_3d,max_c2_3d,n_3d_fine,
5803 3 num_pairs_unlike,num_pairs_unlike_ref)
5804 Call c2_3d(htmp_unlike_3d_coarse,href_unlike_3d_coarse,
5805 1 c2fit_unlike_3d_coarse,c2err_unlike_3d_coarse,
5806 2 max_h_3d,max_c2_3d,n_3d_coarse,
5807 3 num_pairs_unlike,num_pairs_unlike_ref)
5808 else if(mode .eq. 3) then
5809 Call c2_3d(hinc_unlike_3d_fine,href_unlike_3d_fine,
5810 1 c2fit_unlike_3d_fine,c2err_unlike_3d_fine,
5811 2 max_h_3d,max_c2_3d,n_3d_fine,
5812 3 num_pairs_unlike_inc,num_pairs_unlike_ref)
5813 Call c2_3d(hinc_unlike_3d_coarse,href_unlike_3d_coarse,
5814 1 c2fit_unlike_3d_coarse,c2err_unlike_3d_coarse,
5815 2 max_h_3d,max_c2_3d,n_3d_coarse,
5816 3 num_pairs_unlike_inc,num_pairs_unlike_ref)
5819 end if ! End 3D correlations
5825 C-----------------------------------------------------------------------
5828 subroutine c2_1d(h,href,c2,c2err,maxh,maxc2,n,num_pairs_sig,
5832 CCC Computes the two-body correlation function for 1D distributions.
5833 CCC Errors are also computed.
5835 CCC Description of Input Variables in Argument List:
5837 C h(maxh) = signal histogram (numerator)
5838 C href(maxh) = background histogram (denominator)
5839 C c2(maxc2) = correlation function = a/b
5840 C c2err(maxc2) = correlation function error
5841 C maxh = dimension of histogram arrays
5842 C maxc2 = dimension of correlation function array.
5844 C num_pairs_sig = # pairs used in signal histogram
5845 C num_pairs_bkg = # pairs used in background histogram
5848 CCC Local Variable Type Declarations:
5850 integer*4 maxh,maxc2,n,num_pairs_sig,num_pairs_bkg
5851 integer*4 h(maxh), href(maxh)
5854 real*4 c2(maxc2), c2err(maxc2)
5855 real*4 a,a_error,b,b_error
5858 if(href(k).le.0 .or. h(k).le.0) then
5862 a = float(h(k))/float(num_pairs_sig)
5863 a_error = sqrt(float(h(k)))/float(num_pairs_sig)
5864 b = float(href(k))/float(num_pairs_bkg)
5865 b_error = sqrt(float(href(k)))/float(num_pairs_bkg)
5867 c2err(k) = c2(k)*sqrt((a_error/a)**2 + (b_error/b)**2)
5874 C-----------------------------------------------------------------------
5877 subroutine c2_3d(h,href,c2,c2err,maxh,maxc2,n,num_pairs_sig,
5881 CCC Computes the two-body correlation function for 3D distributions.
5882 CCC Errors are also computed.
5884 CCC Description of Input Variables in Argument List:
5886 C h(maxh,maxh,maxh) = 3D signal histogram (numerator)
5887 C href(maxh,maxh,maxh)) = 3D background histogram (denominator)
5888 C c2(maxc2,maxc2,maxc2) = 3D correlation function = a/b
5889 C c2err(maxc2,maxc2,maxc2) = 3D correlation function error
5890 C maxh = dimension of 3D histogram arrays
5891 C maxc2 = dimension of 3D correlation function array.
5893 C num_pairs_sig = # pairs used in signal histogram
5894 C num_pairs_bkg = # pairs used in background histogram
5897 CCC Local Variable Type Declarations:
5899 integer*4 maxh,maxc2,n,num_pairs_sig,num_pairs_bkg
5900 integer*4 h(maxh,maxh,maxh), href(maxh,maxh,maxh)
5903 real*4 c2(maxc2,maxc2,maxc2), c2err(maxc2,maxc2,maxc2)
5904 real*4 a,a_error,b,b_error
5909 if(href(i,j,k).le.0 .or. h(i,j,k).le.0) then
5913 a = float(h(i,j,k))/float(num_pairs_sig)
5914 a_error = sqrt(float(h(i,j,k)))/float(num_pairs_sig)
5915 b = float(href(i,j,k))/float(num_pairs_bkg)
5916 b_error = sqrt(float(href(i,j,k)))/float(num_pairs_bkg)
5918 c2err(i,j,k) = c2(i,j,k)*sqrt((a_error/a)**2 + (b_error/b)**2)
5928 C-------------------------------------------------------------------------
5931 subroutine chisquare(mode,chisq_like_1d,chisq_unlike_1d,
5932 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
5933 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
5934 3 chisq_hist1_1,chisq_hist1_2)
5937 CCC This subroutine calculates the chi-squares for the following:
5938 C o Like pair 1D 2-body correlations
5939 C o Unlike pair 1D 2-body correlations
5940 C o Like pair 3D, Fine Mesh 2-body correlations
5941 C o Unlike pair 3D, Fine Mesh 2-body correlations
5942 C o Like pair 3D, Coarse Mesh 2-body correlations
5943 C o Unlike pair 3D, Coarse Mesh 2-body correlations
5944 C o One-body 1D {pt,phi,eta} (summed) distributions for PID#1
5945 C o One-body 1D {pt,phi,eta} (summed) distributions for PID#2
5947 C (where the separate chi-squares for the 1D pt, phi and eta
5948 C one-body distributions are added and only the sum is returned.)
5950 C 'Mode' determines which one-body histogram is compared to the
5951 C reference histogram, where:
5953 C If mode = 1, then hist1* are used
5954 C If mode = 2, then htmp1* are used
5956 C The one-body reference histograms used in the chi-square calculation
5959 Include 'common_parameters.inc'
5960 Include 'common_mesh.inc'
5961 Include 'common_histograms.inc'
5962 Include 'common_correlations.inc'
5964 CCC Local Variable Type Declarations:
5966 integer*4 i,j,k,mode
5968 real*4 chisq_like_1d, chisq_unlike_1d
5969 real*4 chisq_like_3d_fine,chisq_unlike_3d_fine
5970 real*4 chisq_like_3d_coarse,chisq_unlike_3d_coarse
5971 real*4 chisq_hist1_1,chisq_hist1_2
5973 real*4 n1fac,n2fac ! # part 1(2) used/# part 1(2) used in Ref.
5974 real*4 avgerrsq_pt_1, avgerrsq_phi_1, avgerrsq_eta_1
5975 real*4 avgerrsq_pt_2, avgerrsq_phi_2, avgerrsq_eta_2
5976 real*4 avgerrsq_pt_ref_1,avgerrsq_phi_ref_1,avgerrsq_eta_ref_1
5977 real*4 avgerrsq_pt_ref_2,avgerrsq_phi_ref_2,avgerrsq_eta_ref_2
5980 CCC Initialize all chi-square values to zero:
5983 chisq_unlike_1d = 0.0
5984 chisq_like_3d_fine = 0.0
5985 chisq_unlike_3d_fine = 0.0
5986 chisq_like_3d_coarse = 0.0
5987 chisq_unlike_3d_coarse = 0.0
5991 if(switch_1d .gt. 0) then
5992 if(switch_type.eq.1 .or. switch_type.eq.3) then
5994 if(c2fit_like_1d(i) .ne. 0.0) then
5995 chisq_like_1d = chisq_like_1d + ((c2fit_like_1d(i)
5996 1 - c2mod_like_1d(i))/c2err_like_1d(i))**2
6000 if(switch_type.eq.2 .or. switch_type.eq.3) then
6002 if(c2fit_unlike_1d(i) .ne. 0.0) then
6003 chisq_unlike_1d = chisq_unlike_1d + ((c2fit_unlike_1d(i)
6004 1 - c2mod_unlike_1d(i))/c2err_unlike_1d(i))**2
6008 end if ! End 1D correlation function, chi-square option
6010 if(switch_3d .gt. 0) then
6011 if(switch_type.eq.1 .or. switch_type.eq.3) then
6013 if(n_3d_fine .gt. 0) then
6017 if(c2fit_like_3d_fine(i,j,k).ne.0.0) then
6018 chisq_like_3d_fine = chisq_like_3d_fine
6019 1 + ((c2fit_like_3d_fine(i,j,k)
6020 2 - c2mod_like_3d_fine(i,j,k))
6021 3 /c2err_like_3d_fine(i,j,k))**2
6028 if(n_3d_coarse .gt. 0) then
6029 do i = 1,n_3d_coarse
6030 do j = 1,n_3d_coarse
6031 do k = 1,n_3d_coarse
6032 if((i+j+k).gt.3) then
6033 if(c2fit_like_3d_coarse(i,j,k).ne.0.0) then
6034 chisq_like_3d_coarse = chisq_like_3d_coarse
6035 1 +((c2fit_like_3d_coarse(i,j,k)
6036 2 - c2mod_like_3d_coarse(i,j,k))
6037 3 /c2err_like_3d_coarse(i,j,k))**2
6047 if(switch_type.eq.2 .or. switch_type.eq.3) then
6049 if(n_3d_fine .gt. 0) then
6053 if(c2fit_unlike_3d_fine(i,j,k).ne.0.0) then
6054 chisq_unlike_3d_fine = chisq_unlike_3d_fine
6055 1 + ((c2fit_unlike_3d_fine(i,j,k)
6056 2 - c2mod_unlike_3d_fine(i,j,k))
6057 3 /c2err_unlike_3d_fine(i,j,k))**2
6064 if(n_3d_coarse .gt. 0) then
6065 do i = 1,n_3d_coarse
6066 do j = 1,n_3d_coarse
6067 do k = 1,n_3d_coarse
6068 if((i+j+k).gt.3) then
6069 if(c2fit_unlike_3d_coarse(i,j,k).ne.0.0) then
6070 chisq_unlike_3d_coarse = chisq_unlike_3d_coarse
6071 1 +((c2fit_unlike_3d_coarse(i,j,k)
6072 2 - c2mod_unlike_3d_coarse(i,j,k))
6073 3 /c2err_unlike_3d_coarse(i,j,k))**2
6082 end if ! End of 3D Correlation Function, Chi-Square Option
6084 CCC Obtain chi-squares for one-body distributions
6086 if(pid(1) .gt. 0) then
6087 n1fac = float(n_part_used_1_trk)/float(n_part_used_1_ref)
6088 avgerrsq_pt_1 = float(n_part_used_1_trk)/float(n_pt_bins)
6089 avgerrsq_phi_1 = float(n_part_used_1_trk)/float(n_phi_bins)
6090 avgerrsq_eta_1 = float(n_part_used_1_trk)/float(n_eta_bins)
6091 avgerrsq_pt_ref_1 = float(n_part_used_1_ref)/float(n_pt_bins)
6092 avgerrsq_phi_ref_1 = float(n_part_used_1_ref)/float(n_phi_bins)
6093 avgerrsq_eta_ref_1 = float(n_part_used_1_ref)/float(n_eta_bins)
6096 if(pid(2) .gt. 0) then
6097 n2fac = float(n_part_used_2_trk)/float(n_part_used_2_ref)
6098 avgerrsq_pt_2 = float(n_part_used_2_trk)/float(n_pt_bins)
6099 avgerrsq_phi_2 = float(n_part_used_2_trk)/float(n_phi_bins)
6100 avgerrsq_eta_2 = float(n_part_used_2_trk)/float(n_eta_bins)
6101 avgerrsq_pt_ref_2 = float(n_part_used_2_ref)/float(n_pt_bins)
6102 avgerrsq_phi_ref_2 = float(n_part_used_2_ref)/float(n_phi_bins)
6103 avgerrsq_eta_ref_2 = float(n_part_used_2_ref)/float(n_eta_bins)
6106 if(pid(1) .gt. 0) then
6107 if(mode .eq. 1) then
6110 1 chisq1(hist1_pt_1,href1_pt_1,max_h_1d,avgerrsq_pt_1,
6111 2 avgerrsq_pt_ref_1,n1fac,n_pt_bins)
6112 3 +chisq1(hist1_phi_1,href1_phi_1,max_h_1d,avgerrsq_phi_1,
6113 4 avgerrsq_phi_ref_1,n1fac,n_phi_bins)
6114 5 +chisq1(hist1_eta_1,href1_eta_1,max_h_1d,avgerrsq_eta_1,
6115 6 avgerrsq_eta_ref_1,n1fac,n_eta_bins)
6117 else if(mode .eq. 2) then
6120 1 chisq1(htmp1_pt_1,href1_pt_1,max_h_1d,avgerrsq_pt_1,
6121 2 avgerrsq_pt_ref_1,n1fac,n_pt_bins)
6122 3 +chisq1(htmp1_phi_1,href1_phi_1,max_h_1d,avgerrsq_phi_1,
6123 4 avgerrsq_phi_ref_1,n1fac,n_phi_bins)
6124 5 +chisq1(htmp1_eta_1,href1_eta_1,max_h_1d,avgerrsq_eta_1,
6125 6 avgerrsq_eta_ref_1,n1fac,n_eta_bins)
6128 end if ! End pid(1) one-body histogram chi-square calculation
6130 if(pid(2) .gt. 0) then
6131 if(mode .eq. 1) then
6134 1 chisq1(hist1_pt_2,href1_pt_2,max_h_1d,avgerrsq_pt_2,
6135 2 avgerrsq_pt_ref_2,n2fac,n_pt_bins)
6136 3 +chisq1(hist1_phi_2,href1_phi_2,max_h_1d,avgerrsq_phi_2,
6137 4 avgerrsq_phi_ref_2,n2fac,n_phi_bins)
6138 5 +chisq1(hist1_eta_2,href1_eta_2,max_h_1d,avgerrsq_eta_2,
6139 6 avgerrsq_eta_ref_2,n2fac,n_eta_bins)
6141 else if(mode .eq. 2) then
6144 1 chisq1(htmp1_pt_2,href1_pt_2,max_h_1d,avgerrsq_pt_2,
6145 2 avgerrsq_pt_ref_2,n2fac,n_pt_bins)
6146 3 +chisq1(htmp1_phi_2,href1_phi_2,max_h_1d,avgerrsq_phi_2,
6147 4 avgerrsq_phi_ref_2,n2fac,n_phi_bins)
6148 5 +chisq1(htmp1_eta_2,href1_eta_2,max_h_1d,avgerrsq_eta_2,
6149 6 avgerrsq_eta_ref_2,n2fac,n_eta_bins)
6152 end if ! End pid(2) one-body histogram chi-square calculation
6157 C----------------------------------------------------------------------
6160 real*4 function chisq1(h,href,maxh,herravgsq,hreferravgsq,
6164 CCC Compute chi-square for 1D histogram h(), with respect to the
6165 CCC reference histogram, href().
6167 C h(maxh) = 1D histogram array
6168 C href(maxh) = 1D reference histogram array
6169 C maxh = dimension of histogram arrays
6170 C herravgsq = average error squared in histogram h's bins
6171 C hreferravgsq = average error squared in ref. hist. href's bins
6172 C numfac = ratio of total number of entries in h to that
6174 C nbins = # bins to use in chi-square sum, starting at array
6175 C element 1,2,... nbins (where nbins .le. maxh)
6177 C The chi-square value is returned in chisq1
6179 CCC Local Variable Type Declarations:
6181 integer*4 maxh, nbins, i
6182 integer*4 h(maxh),href(maxh)
6184 real*4 herravgsq,hreferravgsq,numfac,numfacsq
6185 real*4 herrsq,hreferrsq
6188 numfacsq = numfac*numfac
6191 if(h(i) .gt. 0) then
6192 herrsq = float(h(i))
6197 if(href(i) .gt. 0) then
6198 hreferrsq = float(href(i))
6200 hreferrsq = hreferravgsq
6203 chisq1 = chisq1 + ((float(h(i)) - numfac*float(href(i)))**2)
6204 1 /(herrsq + numfacsq*hreferrsq)
6210 C-----------------------------------------------------------------------
6213 Subroutine write_data(mode,ievent)
6216 CCC This subroutine writes the main output file, 'hbt_simulation.out'
6217 C on File Unit 8. File Unit 8 is opened and closed by the main
6220 C Also, the computed 1- and 2-body reference histograms are printed
6221 C out from this subroutine on File Units 11 and 9, respectively. These
6222 C files are opened/closed here.
6224 C Output content determined by input parameter 'mode', where:
6226 C Mode Description of Output
6227 C ----- -----------------------------------------------------------
6228 C 1 basic output file header
6229 C input and derived quantities
6231 C 2 reference histograms (1 and 2-body)
6232 C saved to separate I/O File Unit=11,9 respectively
6234 C 3 reference histogram output
6236 C 4 correlation model
6238 C 5 correlation fit and one-body distributions
6239 C for each event, optional output
6241 C 6 inclusive one-body distributions and inclusive
6242 C correlation fit; projection onto 1D axes.
6245 Include 'common_parameters.inc'
6246 Include 'common_mesh.inc'
6247 Include 'common_histograms.inc'
6248 Include 'common_correlations.inc'
6249 Include 'common_coulomb.inc'
6250 Include 'common_event_summary.inc'
6252 Include 'common_track.inc'
6253 Include 'common_sec_track.inc'
6254 Include 'common_sec_track2.inc'
6256 CCC Local Variable Type Declarations:
6258 integer*4 mode,i,j,k,ievent,ref_print,nev
6260 real*4 nfac1,nfac2,ref_error
6261 real*4 c2mod_proj1(max_c2_3d)
6262 real*4 c2mod_proj2(max_c2_3d)
6263 real*4 c2mod_proj3(max_c2_3d)
6264 real*4 c2fit_proj1(max_c2_3d)
6265 real*4 c2fit_proj2(max_c2_3d)
6266 real*4 c2fit_proj3(max_c2_3d)
6267 real*4 c2err_proj1(max_c2_3d)
6268 real*4 c2err_proj2(max_c2_3d)
6269 real*4 c2err_proj3(max_c2_3d)
6271 C-------------------------------------------
6272 If(mode.eq.1) Then !Basic Output Header
6273 C-------------------------------------------
6278 C write(8,102) n_events
6279 write(8,103) n_pid_types,pid(1),mass1,pid(2),mass2
6280 write(8,104) ref_control
6281 write(8,105) switch_1d
6282 write(8,106) switch_3d
6283 write(8,107) switch_type
6284 write(8,108) switch_coherence
6285 write(8,109) switch_coulomb
6286 write(8,110) switch_fermi_bose
6287 write(8,1101) trk_accep
6288 write(8,111) print_full,print_sector_data
6289 C write(8,112) n_part_used_1_ref,n_part_used_2_ref
6290 C write(8,113) n_part_used_1_inc,n_part_used_2_inc
6291 C write(8,114) num_pairs_like_ref,num_pairs_unlike_ref
6292 C write(8,115) num_pairs_like_inc,num_pairs_unlike_inc
6295 write(8,118) Rside,Rout,Rlong
6296 write(8,119) Rperp,Rparallel,R0
6302 write(8,125) chisq_wt_like_1d
6303 write(8,126) chisq_wt_unlike_1d
6304 write(8,127) chisq_wt_like_3d_fine
6305 write(8,128) chisq_wt_unlike_3d_fine
6306 write(8,129) chisq_wt_like_3d_coarse
6307 write(8,130) chisq_wt_unlike_3d_coarse
6308 write(8,131) chisq_wt_hist1_1
6309 write(8,132) chisq_wt_hist1_2
6311 write(8,134) n_pt_bins,pt_bin_size,pt_min,pt_max
6312 write(8,135) n_phi_bins,phi_bin_size,phi_min,phi_max
6313 write(8,136) n_eta_bins,eta_bin_size,eta_min,eta_max
6315 write(8,138) n_px_bins,delpx,px_min,px_max
6316 write(8,139) n_py_bins,delpy,py_min,py_max
6317 write(8,140) n_pz_bins,delpz,pz_min,pz_max
6318 write(8,141) n_sectors
6320 write(8,143) n_1d_fine,n_1d_coarse,n_1d_total
6321 write(8,144) binsize_1d_fine,binsize_1d_coarse
6322 write(8,145) qmid_1d,qmax_1d
6324 write(8,147) n_3d_fine,n_3d_coarse,n_3d_total
6325 write(8,148) binsize_3d_fine,binsize_3d_coarse
6326 write(8,149) qmid_3d,qmax_3d
6327 write(8,150) n_3d_fine_project
6329 CCC Formats for Mode=1 Output
6331 100 format( 15x,50('*'))
6332 101 format( 15x,'*****',7x,'HBT CORRELATION SIMULATION',7x,'*****')
6333 102 format(///15x,'Number of Events in Event Text Input File=',I5)
6334 103 format( /15x,'#PID types=',I2,' PID#,mass=',I2,F8.5,
6335 1' PID#,mass=',I2,F8.5)
6336 104 format(/ 15x,'Reference Spectra Selection Option=',I2)
6337 105 format(// 15x,'Control Switches: Switch_1d =',I2)
6338 106 format( 15x,' Switch_3d =',I2)
6339 107 format( 15x,' Switch_type =',I2)
6340 108 format( 15x,' Switch_coherence =',I2)
6341 109 format( 15x,' Switch_coulomb =',I2)
6342 110 format( 15x,' Switch_fermi_bose =',I2)
6343 1101 format( 15x,' trk_accep =',F10.7)
6344 111 format(/ 15x,'Print Options: Full=',I2,' Sectors=',I2)
6346 1'Number particles used in Reference, for PID types=',2I5)
6348 1'Number particles used in Inclusive, for PID types=',2I5)
6350 1'Number pairs used in Reference, like and unlike=',2I5)
6352 1'Number pairs used in Inclusive, like and unlike=',2I5)
6353 116 format(// 15x,'Correlation Model Parameters: Chaoticity =',F8.5)
6354 117 format( 15x,'1D Spherical Source Radius=',F8.4)
6355 118 format( 15x,'Bertsch-Pratt R-side,out,long=',3F8.4)
6356 119 format( 15x,'YKP R-perp,parallel,time=',3F8.4)
6357 120 format( 15x,'Coulomb parameter=',F8.5)
6358 121 format(// 15x,'Iteration Controls: Random # seed =',I10)
6359 122 format( 15x,' Max # iterations =',I5)
6360 123 format( 15x,' Momentum Shift Range =',F8.5)
6361 124 format( 15x,' Min % Chi-Sq limit =',F8.5)
6362 125 format(// 15x,'CHI-Sq Weights: correl,like, 1d =',F8.5)
6363 126 format( 15x,' correl,unlike,1d =',F8.5)
6364 127 format( 15x,' correl,like, 3d_fine =',F8.5)
6365 128 format( 15x,' correl,unlike,3d_fine =',F8.5)
6366 129 format( 15x,' correl,like, 3d_coarse=',F8.5)
6367 130 format( 15x,' correl,unlike,3d_coarse=',F8.5)
6368 131 format( 15x,' 1-body, PID#1 =',F8.5)
6369 132 format( 15x,' 1-body, PID#2 =',F8.5)
6370 133 format(// 15x,'Momentum Space Acceptance Range and 1D Bins:')
6371 134 format( 15x,'#bins,bin size,min,max for pt =',I5,3F8.5)
6372 135 format( 15x,'#bins,bin size,min,max for phi=',I5,3F8.4)
6373 136 format( 15x,'#bins,bin size,min,max for eta=',I5,3F8.4)
6374 137 format(// 15x,'Momentum Space Sectors:')
6375 138 format( 15x,'#sectors,sectorsize,min,max for px=',I5,3F8.4)
6376 139 format( 15x,'#sectors,sectorsize,min,max for py=',I5,3F8.4)
6377 140 format( 15x,'#sectors,sectorsize,min,max for pz=',I5,3F8.4)
6378 141 format( 15x,'Total Number of Sectors =',I5)
6379 142 format(// 15x,'2-Body Correlations, 1-D Grid:')
6380 143 format( 15x,'#bins_fine,coarse,total =',3I5)
6381 144 format( 15x,'bin size - fine, coarse =',2F8.5)
6382 145 format( 15x,'Q mid point, Q maximum =',2F8.5)
6383 146 format(// 15x,'2-Body Correlations, 3-D Grid:')
6384 147 format( 15x,'#bins - fine, coarse, total =',3I5)
6385 148 format( 15x,'bin size - fine, coarse =',2F8.5)
6386 149 format( 15x,'Q mid point, Q maximum =',2F8.5)
6387 150 format( 15x,'# 3D fine bin projected =',I5)
6389 CCC END mode=1 Output and Formats
6391 C-----------------------------
6392 Else If(mode.eq.2) Then !Store 2- and 1-body Ref. Histograms
6393 C-----------------------------
6396 open(unit=9,status='unknown',access='sequential',
6397 1 name='hbt_pair_reference.hist')
6398 open(unit=11,status='unknown',access='sequential',
6399 1 name='hbt_singles_reference.hist')
6401 C Write Pair Reference Hist:
6403 write(9,201) n_pid_types,pid(1),pid(2)
6404 write(9,202) n_pt_bins,pt_min,pt_max
6405 write(9,202) n_phi_bins,phi_min,phi_max
6406 write(9,202) n_eta_bins,eta_min,eta_max
6407 write(9,201) switch_1d,switch_3d,switch_type
6408 write(9,203) n_1d_fine,n_1d_coarse,n_3d_fine,n_3d_coarse
6409 write(9,204) binsize_1d_fine,binsize_1d_coarse,
6410 1 binsize_3d_fine,binsize_3d_coarse
6411 write(9,201) num_pairs_like_ref,num_pairs_unlike_ref
6413 202 format(2x,I10,2E15.6)
6415 204 format(2x,4E15.6)
6418 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6419 if(switch_type.eq.1.or.switch_type.eq.3) then
6420 write(9,205) (href_like_1d(i),i=1,n_1d_total)
6422 if(switch_type.eq.2.or.switch_type.eq.3) then
6423 write(9,205) (href_unlike_1d(i),i=1,n_1d_total)
6425 endif !End 1D Ref. Hist. Output
6427 if(switch_3d.gt.0) then
6428 if(switch_type.eq.1.or.switch_type.eq.3) then
6430 if(n_3d_fine.gt.0) then
6434 write(9,205) href_like_3d_fine(i,j,k)
6440 if(n_3d_coarse.gt.0) then
6444 write(9,205) href_like_3d_coarse(i,j,k)
6452 if(switch_type.eq.2.or.switch_type.eq.3) then
6454 if(n_3d_fine.gt.0) then
6458 write(9,205) href_unlike_3d_fine(i,j,k)
6464 if(n_3d_coarse.gt.0) then
6468 write(9,205) href_unlike_3d_coarse(i,j,k)
6475 endif !End 3D Reference Histograms Output
6477 CC Write One-Body - singles histograms:
6479 write(11,201) n_pid_types,pid(1),pid(2)
6480 write(11,202) n_pt_bins,pt_min,pt_max
6481 write(11,202) n_phi_bins,phi_min,phi_max
6482 write(11,202) n_eta_bins,eta_min,eta_max
6483 write(11,201) n_part_used_1_ref,n_part_used_2_ref
6485 if(pid(1).gt.0) then
6486 write(11,205)(href1_pt_1(i),i=1,n_pt_bins)
6487 write(11,205)(href1_phi_1(i),i=1,n_phi_bins)
6488 write(11,205)(href1_eta_1(i),i=1,n_eta_bins)
6492 if(pid(2).gt.0) then
6493 write(11,205)(href1_pt_2(i),i=1,n_pt_bins)
6494 write(11,205)(href1_phi_2(i),i=1,n_phi_bins)
6495 write(11,205)(href1_eta_2(i),i=1,n_eta_bins)
6501 CCC END mode=2 Reference Histogram Output
6503 C-----------------------------
6504 Else If(mode.eq.3) Then !Print out the Reference Histograms
6505 C-----------------------------
6509 write(8,302) n_pt_bins,pt_min,pt_max
6510 write(8,303) n_phi_bins,phi_min,phi_max
6511 write(8,304) n_eta_bins,eta_min,eta_max
6512 write(8,305) n_part_used_1_ref,n_part_used_2_ref
6516 write(8,307) i,href1_pt_1(i),href1_pt_2(i)
6521 write(8,307) i,href1_phi_1(i),href1_phi_2(i)
6526 write(8,307) i,href1_eta_1(i),href1_eta_2(i)
6530 write(8,311) n_1d_fine,n_1d_coarse
6531 write(8,312) binsize_1d_fine,binsize_1d_coarse
6532 write(8,313) n_3d_fine,n_3d_coarse
6533 write(8,314) binsize_3d_fine,binsize_3d_coarse
6534 write(8,315) num_pairs_like_ref,num_pairs_unlike_ref
6536 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6539 write(8,307) i,href_like_1d(i),href_unlike_1d(i)
6541 endif !End Print Out of 2-body, 1D reference histogram
6543 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6548 write(8,318) i,j,k,href_like_3d_fine(i,j,k),
6549 1 href_unlike_3d_fine(i,j,k)
6553 endif !End Print Out of 2-Body, 3D-Fine Mesh Ref. Hist.
6556 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
6561 write(8,318) i,j,k,href_like_3d_coarse(i,j,k),
6562 1 href_unlike_3d_coarse(i,j,k)
6566 endif !End Print Out of 2-Body, 3D-Coarse Mesh Ref. Hist.
6568 CCC Formats for mode=3 Output:
6570 300 format(///5x,15('*'),'REFERENCE HISTOGRAMS',15('*'))
6571 301 format(//15x,'ONE-BODY REFERENCE DISTRIBUTIONS:')
6572 302 format(/ 15x,'PT BINS: (#,min,max)=',I5,2F8.4)
6573 303 format( 15x,'PHI BINS:(#,min,max)=',I5,2F8.4)
6574 304 format( 15x,'ETA BINS:(#,min,max)=',I5,2F8.4)
6575 305 format( 15x,'Number particles used in Ref, PID type1,2=',2I8)
6576 306 format(/ 9x,'PT',10x,'BIN#',5x,'PID-1',5x,'PID-2')
6577 308 format(/ 9x,'PHI',9x,'BIN#',5x,'PID-1',5x,'PID-2')
6578 309 format(/ 9x,'ETA',9x,'BIN#',5x,'PID-1',5x,'PID-2')
6579 307 format( 20x,I5,2I10)
6580 310 format(///15x,'TWO-BODY REFERENCE DISTRIBUTIONS:')
6581 311 format(/ 15x,'#BINS FOR 1D-Fine and Coarse Grid =',2I5)
6582 312 format( 15x,'BIN SIZES FOR 1D-Fine and Coarse =',2F8.5)
6583 313 format( 15x,'#BINS FOR 3D-Fine and Coarse Grid =',2I5)
6584 314 format( 15x,'BIN SIZES FOR 3D-Fine and Coarse =',2F8.5)
6585 315 format( 15x,'Number of Like and Unlike Pairs For Ref. = ',
6587 316 format(/5x,'2-BODY, 1D',6x,'BIN#',5x,'LIKE',5x,'UNLIKE')
6588 317 format(/3x,'2-BODY, 3D-FINE',2x,'BIN:i',4x,'j',4x,
6589 1 'k',5x,'LIKE',5x,'UNLIKE')
6590 318 format( 20x,3I5,2I10)
6591 319 format(/2x,'2-BODY, 3D-COARSE',1x,'BIN:i',4x,'j',4x,
6592 1 'k',5x,'LIKE',5x,'UNLIKE')
6594 CCC END mode=3 Output and Formats
6596 C----------------------------
6597 Else If(mode.eq.4) Then !Print Correlation Function Model
6598 C----------------------------
6601 write(8,311) n_1d_fine,n_1d_coarse
6602 write(8,312) binsize_1d_fine,binsize_1d_coarse
6603 write(8,313) n_3d_fine,n_3d_coarse
6604 write(8,314) binsize_3d_fine,binsize_3d_coarse
6607 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6610 write(8,407) i,c2mod_like_1d(i),c2mod_unlike_1d(i)
6612 endif !End Print Out of 2-body, 1D Model Correction Functions
6614 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6619 write(8,418) i,j,k,c2mod_like_3d_fine(i,j,k),
6620 1 c2mod_unlike_3d_fine(i,j,k)
6624 endif !End Print Out of 2-Body, 3D-Fine mesh Model Correl. Function
6626 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
6631 write(8,418) i,j,k,c2mod_like_3d_coarse(i,j,k),
6632 1 c2mod_unlike_3d_coarse(i,j,k)
6636 endif !End Print Out of 2-Body, 3D-Coarse Model Correlation Function
6638 if(switch_coulomb.eq.3) then ! Print interpolated Pratt Model
6639 CC ! Coulomb Correction for Finite
6640 CC ! Source Radius Q0
6644 write(8,403) i,q_coul(i),c2_coul_like(i),
6649 CCC Additional Formats for Mode=4 Output:
6651 400 format(///5x,15('*'),'MODEL CORRELATION FUNCTIONS',15('*'))
6652 401 format(///15x,'COULOMB SOURCE RADIUS FOR PRATT MODEL=',F8.4)
6653 402 format(// 5x,'q-bin',2x,'q',4x,'C2_coul_like',2x,
6655 403 format( 5x,I5,3E15.6)
6656 407 format( 20x,I5,2F10.7)
6657 418 format( 20x,3I5,2F10.7)
6659 CCC END MODE = 4 OUTPUT AND FORMATS
6661 C------------------------------
6662 Else If(mode.eq.5) Then ! Optional Output for 1- and 2-Body Fits
6663 C------------------------------ ! for each event.
6666 write(8,501) n_part_1_trk,n_part_2_trk,n_part_tot_trk
6667 write(8,502) n_part_used_1_trk, n_part_used_2_trk
6668 write(8,503) num_pairs_like, num_pairs_unlike
6670 CCC Output one-body distributions for event:
6673 if(pid(1) .gt. 0) then
6674 nfac1 = float(n_part_used_1_trk)/float(n_part_used_1_ref)
6679 ref_print = int(nfac1*float(href1_pt_1(i)))
6680 ref_error = nfac1*sqrt(float(href1_pt_1(i)))
6681 write(8,510) i,hist1_pt_1(i),ref_print,ref_error
6686 ref_print = int(nfac1*float(href1_phi_1(i)))
6687 ref_error = nfac1*sqrt(float(href1_phi_1(i)))
6688 write(8,510) i,hist1_phi_1(i),ref_print,ref_error
6693 ref_print = int(nfac1*float(href1_eta_1(i)))
6694 ref_error = nfac1*sqrt(float(href1_eta_1(i)))
6695 write(8,510) i,hist1_eta_1(i),ref_print,ref_error
6698 end if ! End PID # 1, One-Body Distribution Output
6700 if(pid(2) .gt. 0) then
6701 nfac2 = float(n_part_used_2_trk)/float(n_part_used_2_ref)
6706 ref_print = int(nfac2*float(href1_pt_2(i)))
6707 ref_error = nfac2*sqrt(float(href1_pt_2(i)))
6708 write(8,510) i,hist1_pt_2(i),ref_print,ref_error
6713 ref_print = int(nfac2*float(href1_phi_2(i)))
6714 ref_error = nfac2*sqrt(float(href1_phi_2(i)))
6715 write(8,510) i,hist1_phi_2(i),ref_print,ref_error
6720 ref_print = int(nfac2*float(href1_eta_2(i)))
6721 ref_error = nfac2*sqrt(float(href1_eta_2(i)))
6722 write(8,510) i,hist1_eta_2(i),ref_print,ref_error
6725 end if ! End PID # 2, One-Body Distribution Output
6727 CCC Output Two-Body Correlation Functions for Event:
6730 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
6735 write(8,523) i,c2mod_like_1d(i),c2fit_like_1d(i),
6736 1 c2err_like_1d(i),c2mod_unlike_1d(i),
6737 2 c2fit_unlike_1d(i),c2err_unlike_1d(i)
6739 end if ! End 1D Correlation Model and Fit Output
6741 if(switch_3d.gt.0 .and. n_3d_fine.gt.0) then
6748 write(8,526) i,j,k,c2mod_like_3d_fine(i,j,k),
6749 1 c2fit_like_3d_fine(i,j,k),c2err_like_3d_fine(i,j,k),
6750 2 c2mod_unlike_3d_fine(i,j,k),c2fit_unlike_3d_fine(i,j,k),
6751 3 c2err_unlike_3d_fine(i,j,k)
6755 end if ! End 3D Fine Mesh Correlation Model and Fit Output
6757 if(switch_3d.gt.0 .and. n_3d_coarse.gt.0) then
6761 do i = 1,n_3d_coarse
6762 do j = 1,n_3d_coarse
6763 do k = 1,n_3d_coarse
6764 write(8,526) i,j,k,c2mod_like_3d_coarse(i,j,k),
6765 1 c2fit_like_3d_coarse(i,j,k),c2err_like_3d_coarse(i,j,k),
6766 2 c2mod_unlike_3d_coarse(i,j,k),c2fit_unlike_3d_coarse(i,j,k),
6767 3 c2err_unlike_3d_coarse(i,j,k)
6771 end if ! End 3D Coarse Mesh Correlation Model and Fit Output
6773 CCC Output Event Summary and Chi-Square Information for Event:
6776 write(8,540) num_iter(ievent)
6777 write(8,541) n_part_used_1_store(ievent),
6778 1 n_part_used_2_store(ievent)
6779 write(8,5411) n_part_tot_store(ievent)
6780 write(8,542) num_sec_flagged_store(ievent)
6781 write(8,543) frac_trks_out(ievent),frac_trks_flag(ievent)
6782 write(8,544) chisq_like_1d_store(ievent),
6783 1 chisq_unlike_1d_store(ievent)
6784 write(8,545) chisq_like_3d_fine_store(ievent),
6785 1 chisq_unlike_3d_fine_store(ievent)
6786 write(8,546) chisq_like_3d_coarse_store(ievent),
6787 1 chisq_unlike_3d_coarse_store(ievent)
6788 write(8,547) chisq_hist1_1_store(ievent),
6789 1 chisq_hist1_2_store(ievent)
6790 write(8,548) chisq_total_store(ievent)
6792 CCC Formats for Mode = 5 Output:
6794 500 Format(///5x,5('*'),'Fitted 1-Body Distributions and ',
6795 1 'Correlations for Event #',I5,5('*'))
6796 501 Format(//15x,'Number of Particles of PID types 1,2,total = ',
6798 502 Format( 15x,'Number of Particles of PID types 1,2 Used = ',
6800 503 Format( 15x,'Number of Pairs Used - Like and Unlike = ',2I10)
6801 504 Format(//5x,'Fitted and Normalized Reference One-Body ',
6803 505 Format( /10x,'Particle Type 1: Reference Scale Factor = ',E12.5)
6804 506 Format( /10x,'Particle Type 2: Reference Scale Factor = ',E12.5)
6805 507 Format(/2x,' PT: BIN#',5x,'hist1',7x,'href1-scaled',2x,
6807 508 Format(/2x,'PHI: BIN#',5x,'hist1',7x,'href1-scaled',2x,
6809 509 Format(/2x,'ETA: BIN#',5x,'hist1',7x,'href1-scaled',2x,
6811 510 Format(7x,I4,3x,I7,8x,I7,7x,F10.5)
6812 520 Format(//5x,'Model and Fitted Correlations')
6813 530 Format(//21x,'One-Dimensional Fit - Fine & Coarse Mesh')
6814 531 Format(//25x,'Three-Dimensional Fit - Fine Mesh')
6815 532 Format(//24x,'Three-Dimensional Fit - Coarse Mesh')
6816 521 Format(/1x,'BIN',13x,'LIKE PAIRS',27x,'UNLIKE PAIRS')
6817 522 Format(8x,'MOD',9x,'FIT',9x,'ERR',11x,'MOD',9x,'FIT',9x,
6819 523 Format(1x,I3,3E12.4,2x,3E12.4)
6820 524 Format(/2x,'BINS',12x,'LIKE PAIRS',24x,'UNLIKE PAIRS')
6821 525 Format(1x,' i j k',4x,'MOD',8x,'FIT',8x,'ERR',10x,'MOD',
6822 1 8x,'FIT',8x,'ERR',/)
6823 526 Format(1x,3I2,3E11.4,2x,3E11.4)
6824 539 Format(///10x,'Event and Chi-Square Summary for Event #',I5)
6825 540 Format( //15x,'Number of Iterations =',F10.2)
6826 541 Format( 15x,'# Particles Used for PID Types1,2=',2F10.2)
6827 5411 Format( 15x,'Total # Particles in track table =',F10.2)
6828 542 Format( 15x,'# Sectors Flagged =',F10.2)
6829 543 Format( 15x,'Frac Trks Out of Accep., Flagged =',2E11.4)
6830 544 Format( 15x,'Chi-Sq: 1D - Like & Unlike =',2E11.4)
6831 545 Format( 15x,'Chi-Sq: 3D - Fine -Like & Unlike =',2E11.4)
6832 546 Format( 15x,'Chi-Sq: 3D - Coarse-Like &Unlike =',2E11.4)
6833 547 Format( 15x,'Chi-Sq: One-Body Dist. PID# 1&2 =',2E11.4)
6834 548 Format( 15x,'Chi-Sq: Total Weighted =',E11.4)
6836 CCC End Mode = 5 Output and Formats
6838 C------------------------------
6839 Else If(mode.eq.6) Then ! Inclusive 1 & 2 Body Output
6840 C------------------------------
6841 write(8,600) n_events
6842 write(8,601) n_part_used_1_inc,n_part_used_2_inc
6843 write(8,602) num_pairs_like_inc,num_pairs_unlike_inc
6846 if(pid(1).gt.0) then
6847 C Division by zero check
6848 IF (n_part_used_1_ref .LE. 0) THEN
6849 PRINT*,'************************************'
6850 PRINT*,'* HBT PROCESSOR *'
6851 PRINT*,'* Number of particles selected for *'
6852 PRINT*,'* processing is less or equal *'
6853 PRINT*,'* !!!!!!!! ZER0 !!!!!!!!!! *'
6854 PRINT*,'* unable to proceed *'
6855 PRINT*,'* EXITING FORTRAN *'
6857 PRINT*,'* HINT: broad the parameter regions*'
6858 PRINT*,'* OR/AND number of particles OR/AND*'
6859 PRINT*,'* number of events *'
6860 PRINT*,'************************************'
6862 5481 FORMAT(5x,'Number of particles selected for processing is',
6863 1 ' less or equal 0',
6868 nfac1=float(n_part_used_1_inc)/float(n_part_used_1_ref)
6872 ref_print=int(nfac1*float(href1_pt_1(i)))
6873 ref_error=nfac1*sqrt(float(href1_pt_1(i)))
6874 write(8,510) i,hinc1_pt_1(i),ref_print,ref_error
6879 ref_print=int(nfac1*float(href1_phi_1(i)))
6880 ref_error=nfac1*sqrt(float(href1_phi_1(i)))
6881 write(8,510) i,hinc1_phi_1(i),ref_print,ref_error
6886 ref_print=int(nfac1*float(href1_eta_1(i)))
6887 ref_error=nfac1*sqrt(float(href1_eta_1(i)))
6888 write(8,510) i,hinc1_eta_1(i),ref_print,ref_error
6891 endif !END PID #1 One-BODY INCL. DISTRIBUTION OUTPUT
6893 if(pid(2).gt.0) then
6894 nfac2=float(n_part_used_2_inc)/float(n_part_used_2_ref)
6898 ref_print=int(nfac2*float(href1_pt_2(i)))
6899 ref_error=nfac2*sqrt(float(href1_pt_2(i)))
6900 write(8,510) i,hinc1_pt_2(i),ref_print,ref_error
6905 ref_print=int(nfac2*float(href1_phi_2(i)))
6906 ref_error=nfac2*sqrt(float(href1_phi_2(i)))
6907 write(8,510) i,hinc1_phi_2(i),ref_print,ref_error
6912 ref_print=int(nfac2*float(href1_eta_2(i)))
6913 ref_error=nfac2*sqrt(float(href1_eta_2(i)))
6914 write(8,510) i,hinc1_eta_2(i),ref_print,ref_error
6917 endif !END PID #2 One-BODY INCl. DISTRIBUTION OUTPUT
6919 CC OUTPUT TWO-BODY INCLUSIVE HISTOGRAMS:
6922 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6925 write(8,307) i,hinc_like_1d(i),hinc_unlike_1d(i)
6927 end if ! End Print out of 2-Body, 1D Inclusive Histograms
6929 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6934 write(8,318) i,j,k,hinc_like_3d_fine(i,j,k),
6935 1 hinc_unlike_3d_fine(i,j,k)
6939 endif ! End Print out of 2-Body, 3D-Fine Inclusive Histograms
6941 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
6946 write(8,318) i,j,k,hinc_like_3d_coarse(i,j,k),
6947 1 hinc_unlike_3d_coarse(i,j,k)
6951 endif ! End Print out of 2-Body, 3D-Coarse Inclusive Histograms
6953 CC OUTPUT TWO-BODY INCL.CORRELATION FUNCTIONS FOR EVENT
6956 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6961 write(8,523) i,c2mod_like_1d(i),c2fit_like_1d(i),
6962 1 c2err_like_1d(i),c2mod_unlike_1d(i),
6963 2 c2fit_unlike_1d(i),c2err_unlike_1d(i)
6967 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6974 write(8,526) i,j,k,c2mod_like_3d_fine(i,j,k),
6975 1 c2fit_like_3d_fine(i,j,k),
6976 2 c2err_like_3d_fine(i,j,k),
6977 3 c2mod_unlike_3d_fine(i,j,k),
6978 4 c2fit_unlike_3d_fine(i,j,k),
6979 5 c2err_unlike_3d_fine(i,j,k)
6986 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
6993 write(8,526) i,j,k,c2mod_like_3d_coarse(i,j,k),
6994 1 c2fit_like_3d_coarse(i,j,k),
6995 2 c2err_like_3d_coarse(i,j,k),
6996 3 c2mod_unlike_3d_coarse(i,j,k),
6997 4 c2fit_unlike_3d_coarse(i,j,k),
6998 5 c2err_unlike_3d_coarse(i,j,k)
7005 CCC Compute and Print 1D projections of 3D fine mesh C2 model,
7006 CCC fit and errors for like and unlike pairs.
7008 if(switch_3d .gt. 0 .and. n_3d_fine .gt. 0) then
7009 if(switch_type .eq. 1 .or. switch_type .eq. 3) then
7010 Call c2_3d_projected(hinc_like_3d_fine,
7011 1 href_like_3d_fine,c2mod_like_3d_fine,
7012 2 c2mod_proj1,c2mod_proj2,c2mod_proj3,
7013 3 c2fit_proj1,c2fit_proj2,c2fit_proj3,
7014 4 c2err_proj1,c2err_proj2,c2err_proj3,
7015 5 max_h_3d, max_c2_3d, n_3d_fine,
7016 6 n_3d_fine_project,num_pairs_like_inc,
7017 7 num_pairs_like_ref)
7022 write(8,658) i,c2mod_proj1(i),c2fit_proj1(i),c2err_proj1(i)
7027 write(8,658) i,c2mod_proj2(i),c2fit_proj2(i),c2err_proj2(i)
7032 write(8,658) i,c2mod_proj3(i),c2fit_proj3(i),c2err_proj3(i)
7034 end if ! End Like pair output
7036 if(switch_type .eq. 2 .or. switch_type .eq. 3) then
7037 Call c2_3d_projected(hinc_unlike_3d_fine,
7038 1 href_unlike_3d_fine,c2mod_unlike_3d_fine,
7039 2 c2mod_proj1,c2mod_proj2,c2mod_proj3,
7040 3 c2fit_proj1,c2fit_proj2,c2fit_proj3,
7041 4 c2err_proj1,c2err_proj2,c2err_proj3,
7042 5 max_h_3d, max_c2_3d, n_3d_fine,
7043 6 n_3d_fine_project,num_pairs_unlike_inc,
7044 7 num_pairs_unlike_ref)
7048 write(8,658) i,c2mod_proj1(i),c2fit_proj1(i),c2err_proj1(i)
7053 write(8,658) i,c2mod_proj2(i),c2fit_proj2(i),c2err_proj2(i)
7058 write(8,658) i,c2mod_proj3(i),c2fit_proj3(i),c2err_proj3(i)
7060 end if ! End Unlike pair output
7061 end if ! End 1D projections
7064 CCC EVENT AND CHISQ SUMMARY INFORMATION:
7066 if(n_events.le.max_events) then
7076 write(8,623) i,num_iter(i),n_part_used_1_store(i),
7077 1 n_part_used_2_store(i),
7078 2 num_sec_flagged_store(i),
7079 3 frac_trks_out(i),frac_trks_flag(i),
7080 4 chisq_total_store(i)
7083 write(8,6231) trk_maxlen
7086 write(8,6233) i,n_part_tot_store(i)
7091 write(8,625) i,chisq_like_1d_store(i),
7092 1 chisq_unlike_1d_store(i)
7098 write(8,625) i,chisq_like_3d_fine_store(i),
7099 1 chisq_unlike_3d_fine_store(i)
7105 write(8,625) i,chisq_like_3d_coarse_store(i),
7106 1 chisq_unlike_3d_coarse_store(i)
7112 write(8,625) i,chisq_hist1_1_store(i),
7113 1 chisq_hist1_2_store(i)
7116 CCC Output the Mean and RMS values for the Event Loop:
7120 write(8,631) niter_mean,niter_rms
7121 write(8,632) npart1_mean,npart1_rms
7122 write(8,633) npart2_mean,npart2_rms
7123 write(8,6331) npart_tot_mean, npart_tot_rms
7124 write(8,634) nsec_flag_mean,nsec_flag_rms
7125 write(8,635) frac_trks_out_mean,frac_trks_out_rms
7126 write(8,636) frac_trks_flag_mean,frac_trks_flag_rms
7127 write(8,637) chi_l1d_mean,chi_l1d_rms
7128 write(8,638) chi_u1d_mean,chi_u1d_rms
7129 write(8,639) chi_l3f_mean,chi_l3f_rms
7130 write(8,640) chi_u3f_mean,chi_u3f_rms
7131 write(8,641) chi_l3c_mean,chi_l3c_rms
7132 write(8,642) chi_u3c_mean,chi_u3c_rms
7133 write(8,643) chi_1_1_mean,chi_1_1_rms
7134 write(8,644) chi_1_2_mean,chi_1_2_rms
7135 write(8,645) chi_tot_mean,chi_tot_rms
7137 CCC FORMATS FOR MODE = 6 OUTPUT
7139 600 format(/// 2x,'FITTED 1-BODY DIST. AND CORRELATIONS ',
7140 1 'FOR INCLUSIVE SUM OF',I5,' EVENTS')
7141 601 format(// 15x,'Inclusive # Particles USED of PID ',
7143 602 format( 15x,'Inclusive # of pairs used; like/unlike=',
7145 603 format(// 5x,'Inclusive and Normalized Reference ',
7146 1 'One-Body Distributions')
7147 604 format(/ 10x,'Inclusive: Particle Type 1 - Reference ',
7148 1 'Scale Factor=',E12.5)
7149 605 format(/ 2x,' PT: BIN#',5x,'hinc1',7x,'href1-scaled',2x,
7151 606 format(/ 2x,'PHI: BIN#',5x,'hinc1',7x,'href1-scaled',2x,
7153 607 format(/ 2x,'ETA: BIN#',5x,'hinc1',7x,'href1-scaled',2x,
7155 608 format(/ 10x,'Inclusive: Particle Type 2 - ',
7156 1 'Reference Scale Factor=',E12.5)
7157 620 format(// 5x,'MODEL AND INCLUSIVE FITTED CORRELATIONS')
7158 621 format(// 15x,'Event and Chi-Square Summary Lists')
7159 622 format(/ 3x,'event',2x,'#iter',3x,'#PID1',4x,'#PID2',3x,
7160 1 '#sec-flg',3x,'frac-out',4x,'frac-flg',3x,'CHISQ-TOT')
7161 623 format(3x,I5,2x,F6.0,2(1x,F8.0),1x,F9.0,
7162 1 2(1x,F11.8),1x,E11.4)
7163 6231 format(/5x,'Max# tracks allowed in track table = ',I8)
7164 6232 format(/5x,'event',4x,'Tot# trks')
7165 6233 format(5x,I5,F12.2)
7166 624 format(/5x,'event',4x,'CHI_l1d',8x,'CHI_u1d')
7167 626 format(/5x,'event',4x,'CHI_l3f',8x,'CHI_u3f')
7168 627 format(/5x,'event',4x,'CHI_l3c',8x,'CHI_u3c')
7169 628 format(/5x,'event',4x,'CHI_1-1',8x,'CHI_1-2')
7170 625 format(5x,I5,2E15.6)
7171 629 format(// 10x,'Event and Chi-Square Summary - ',
7172 1 'Mean and RMS Values')
7173 630 format(/ 14x,'Quantity',15x,'Mean',11x,'RMS')
7174 631 format( 5x,'Number of Iterations ',2E15.6)
7175 632 format( 5x,'#PID Type 1 ',2E15.6)
7176 633 format( 5x,'#PID Type 2 ',2E15.6)
7177 6331 format( 5x,'Tot # Tracks in Table ',2E15.6)
7178 634 format( 5x,'#Sectors Flagged ',2E15.6)
7179 635 format( 5x,'Frac. Trks Out of Accept. ',2E15.6)
7180 636 format( 5x,'Frac. Trks Flagged ',2E15.6)
7181 637 format( 5x,'CHISQ like 1D ',2E15.6)
7182 638 format( 5x,'CHISQ unlike 1D ',2E15.6)
7183 639 format( 5x,'CHISQ like 3D Fine ',2E15.6)
7184 640 format( 5x,'CHISQ unlike 3D Fine ',2E15.6)
7185 641 format( 5x,'CHISQ like 3D Coarse ',2E15.6)
7186 642 format( 5x,'CHISQ unlike 3D Coarse ',2E15.6)
7187 643 format( 5x,'CHISQ 1 Body #1 ',2E15.6)
7188 644 format( 5x,'CHISQ 1 Body #2 ',2E15.6)
7189 645 format( 5x,'CHISQ Total ',2E15.6)
7190 650 format(//10x ,'Inclusive Three-Dimensional Projected Fits -',
7192 651 format( /25x ,'Like Pairs - Axis #1 ')
7193 652 format( /25x ,'Like Pairs - Axis #2 ')
7194 653 format( /25x ,'Like Pairs - Axis #3 ')
7195 654 format( /25x ,'Unlike Pairs - Axis #1 ')
7196 655 format( /25x ,'Unlike Pairs - Axis #2 ')
7197 656 format( /25x ,'Unlike Pairs - Axis #3 ')
7198 657 format( 2x,'BIN#',3x,'Model',8x,'Fit',8x,'Error')
7199 658 format( 3x,I3,3E12.4)
7200 660 format(// 5x,'INCLUSIVE TWO-BODY HISTOGRAMS')
7202 CCC END MODE = 6 OUTPUT AND FORMATS
7211 C-----------------------------------------------------------------------
7214 subroutine c2_3d_projected(h,href,c2mod,
7215 1 c2mod_proj1,c2mod_proj2,c2mod_proj3,
7216 2 c2fit_proj1,c2fit_proj2,c2fit_proj3,
7217 3 c2err_proj1,c2err_proj2,c2err_proj3,
7218 4 maxh,maxc2,n,n_proj,num_pairs_sig,
7223 CCC This Subroutine computes the projected two-body correlation
7224 CCC function for 3D distributions - fine mesh only; for both the
7225 CCC correlation model (weighted with the reference histogram) and
7226 CCC the inclusive correlation fit.
7228 CCC Description of Input Variables in the Argument List:
7230 CCC h(maxh,maxh,maxh) = 3D fine mesh inclusive signal histog.
7231 CCC href(maxh,maxh,maxh) = 3D fine mesh inclusive background hist.
7232 CCC c2mod(maxc2,maxc2,maxc2) = 3D fine mesh correlation model
7233 CCC maxh = Dimension of 3D fine mesh histogram arrays
7234 CCC maxc2 = Dimension of 3D fine mesh correlation function arrays
7235 CCC n = Number of bins to use
7236 CCC n_proj = Number of bins to integrate in (i,j) to project onto (k)
7237 CCC num_pairs_sig = # pairs used in signal histogram
7238 CCC num_pairs_bkg = # pairs used in background histogram
7240 CCC Description of Output quantities:
7242 CCC c2mod_proj1,2,3(maxc2) = Reference histogram weighted 1D projections
7243 CCC of C2 model function along {1,2,3} axes.
7244 CCC c2fit_proj1,2,3(maxc2) = Fitted 3D correlation function projected
7245 CCC onto {1,2,3} axes.
7246 CCC c2err_proj1,2,3(maxc2) = Error in fitted 3D correlation function
7247 CCC projected onto {1,2,3} axes.
7249 CCC Local Variable Type Declarations:
7251 integer*4 maxh,maxc2,n,n_proj,num_pairs_sig,num_pairs_bkg
7252 integer*4 h(maxh,maxh,maxh),href(maxh,maxh,maxh)
7255 real*4 c2mod(maxc2,maxc2,maxc2)
7256 real*4 c2mod_proj1(maxc2),c2mod_proj2(maxc2),c2mod_proj3(maxc2)
7257 real*4 c2fit_proj1(maxc2),c2fit_proj2(maxc2),c2fit_proj3(maxc2)
7258 real*4 c2err_proj1(maxc2),c2err_proj2(maxc2),c2err_proj3(maxc2)
7259 real*4 a,a_error,b,b_error
7260 real*4 sum1n,sum1d,sum2n,sum2d,sum3n,sum3d
7262 CCC Initialize arrays to zero:
7265 c2mod_proj1(i) = 0.0
7266 c2mod_proj2(i) = 0.0
7267 c2mod_proj3(i) = 0.0
7268 c2fit_proj1(i) = 0.0
7269 c2fit_proj2(i) = 0.0
7270 c2fit_proj3(i) = 0.0
7271 c2err_proj1(i) = 0.0
7272 c2err_proj2(i) = 0.0
7273 c2err_proj3(i) = 0.0
7276 CCC Project Reference spectra (histogram) weighted model correlation:
7287 sum1n = sum1n + c2mod(i,j,k)*float(href(i,j,k))
7288 sum1d = sum1d + float(href(i,j,k))
7289 sum2n = sum2n + c2mod(k,i,j)*float(href(k,i,j))
7290 sum2d = sum2d + float(href(k,i,j))
7291 sum3n = sum3n + c2mod(j,k,i)*float(href(j,k,i))
7292 sum3d = sum3d + float(href(j,k,i))
7295 if(sum1d .le. 0.0) then
7296 c2mod_proj1(i) = 0.0
7298 c2mod_proj1(i) = sum1n/sum1d
7300 if(sum2d .le. 0.0) then
7301 c2mod_proj2(i) = 0.0
7303 c2mod_proj2(i) = sum2n/sum2d
7305 if(sum3d .le. 0.0) then
7306 c2mod_proj3(i) = 0.0
7308 c2mod_proj3(i) = sum3n/sum3d
7312 CCC Calculate and Project the fitted correlation functions:
7323 sum1n = sum1n + float(h(i,j,k))
7324 sum1d = sum1d + float(href(i,j,k))
7325 sum2n = sum2n + float(h(k,i,j))
7326 sum2d = sum2d + float(href(k,i,j))
7327 sum3n = sum3n + float(h(j,k,i))
7328 sum3d = sum3d + float(href(j,k,i))
7331 if(sum1n .le. 0.0 .or. sum1d .le. 0.0) then
7332 c2fit_proj1(i) = 0.0
7333 c2err_proj1(i) = 1.0
7335 a = sum1n/float(num_pairs_sig)
7336 a_error = sqrt(sum1n)/float(num_pairs_sig)
7337 b = sum1d/float(num_pairs_bkg)
7338 b_error = sqrt(sum1d)/float(num_pairs_bkg)
7339 c2fit_proj1(i) = a/b
7340 c2err_proj1(i) = c2fit_proj1(i)*sqrt((a_error/a)**2
7343 if(sum2n .le. 0.0 .or. sum2d .le. 0.0) then
7344 c2fit_proj2(i) = 0.0
7345 c2err_proj2(i) = 1.0
7347 a = sum2n/float(num_pairs_sig)
7348 a_error = sqrt(sum2n)/float(num_pairs_sig)
7349 b = sum2d/float(num_pairs_bkg)
7350 b_error = sqrt(sum2d)/float(num_pairs_bkg)
7351 c2fit_proj2(i) = a/b
7352 c2err_proj2(i) = c2fit_proj2(i)*sqrt((a_error/a)**2
7355 if(sum3n .le. 0.0 .or. sum3d .le. 0.0) then
7356 c2fit_proj3(i) = 0.0
7357 c2err_proj3(i) = 1.0
7359 a = sum3n/float(num_pairs_sig)
7360 a_error = sqrt(sum3n)/float(num_pairs_sig)
7361 b = sum3d/float(num_pairs_bkg)
7362 b_error = sqrt(sum3d)/float(num_pairs_bkg)
7363 c2fit_proj3(i) = a/b
7364 c2err_proj3(i) = c2fit_proj3(i)*sqrt((a_error/a)**2
7372 C----------------------------------------------------------------------
7374 C>>>>>>>>>>>>>> Piotr, this needs to be replaced with Ali random
7375 C>>>>>>>>>>>>>>> number generator
7377 * real*4 function hbtpran(i)
7381 * Call ranhbtp(r,1,i)
7386 * Include 'ranlux2.f'
7387 C----------------------------------------------------------------------