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:
2720 C write(*,*) 'REF HISTO N Ev = ', n_events
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)
2735 if (href1_pt_1(ipt) .lt. 0 ) then
2736 write(*,*) 'href1_pt_1 bin ',ipt,'is less then 0'
2740 do iphi = 1,n_phi_bins
2741 href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2742 if (href1_phi_1(iphi) .lt. 0 ) then
2743 write(*,*) 'href1_phi_1 bin ',iphi,'is less then 0'
2747 do ieta = 1,n_eta_bins
2748 href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2749 if (href1_eta_1(ieta) .lt. 0 ) then
2750 write(*,*) 'href1_eta_1 bin ',ieta,'is less then 0'
2755 if(pid(2) .gt. 0) then
2756 Call histog1(1,0,2,pid(2),0.0,0.0,0.0)
2757 n_part_used_2_ref = n_part_used_2_ref + n_part_used_2_trk
2759 do ipt = 1,n_pt_bins
2760 href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2763 do iphi = 1,n_phi_bins
2764 href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2767 do ieta = 1,n_eta_bins
2768 href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2772 else if(sign_toggle .eq. (-1)) then ! Put tracks into 'trk2'
2775 Call stm_build(2,0,0)
2776 if(pid(1) .gt. 0) then
2777 Call histog1(4,0,1,pid(1),0.0,0.0,0.0)
2778 n_part_used_1_ref = n_part_used_1_ref +n_part_used_1_trk2
2780 do ipt = 1,n_pt_bins
2781 href1_pt_1(ipt) = href1_pt_1(ipt) + hist1_pt_1(ipt)
2784 do iphi = 1,n_phi_bins
2785 href1_phi_1(iphi) = href1_phi_1(iphi) + hist1_phi_1(iphi)
2788 do ieta = 1,n_eta_bins
2789 href1_eta_1(ieta) = href1_eta_1(ieta) + hist1_eta_1(ieta)
2793 if(pid(2) .gt. 0) then
2794 Call histog1(4,0,2,pid(2),0.0,0.0,0.0)
2795 n_part_used_2_ref = n_part_used_2_ref +n_part_used_2_trk2
2797 do ipt = 1,n_pt_bins
2798 href1_pt_2(ipt) = href1_pt_2(ipt) + hist1_pt_2(ipt)
2801 do iphi = 1,n_phi_bins
2802 href1_phi_2(iphi) = href1_phi_2(iphi) + hist1_phi_2(iphi)
2805 do ieta = 1,n_eta_bins
2806 href1_eta_2(ieta) = href1_eta_2(ieta) + hist1_eta_2(ieta)
2810 end if ! End read and load to trk or trk2 option
2812 sign_toggle = -sign_toggle
2814 if(i .gt. 1) then ! Compute 2-body reference histograms
2815 Call histog2(4,0,0,0,0,0.0,0.0,0.0,0.0)
2816 num_pairs_like_ref = num_pairs_like_ref
2817 1 + n_part_used_1_trk * n_part_used_1_trk2
2818 2 + n_part_used_2_trk * n_part_used_2_trk2
2819 num_pairs_unlike_ref = num_pairs_unlike_ref
2820 1 + n_part_used_1_trk * n_part_used_2_trk2
2821 2 + n_part_used_2_trk * n_part_used_1_trk2
2822 C write(*,*) 'num_pairs_like_ref',num_pairs_like_ref
2823 C write(*,*) 'num_pairs_unlike_ref',num_pairs_unlike_ref
2826 end do ! End of Event Loop
2828 CCC Write out the pair and one-body reference Histograms:
2829 Call write_data(2,0)
2831 If(ALICE .eq. 0) Then
2836 end if ! End Reference Histogram read/calculate option
2841 C----------------------------------------------------------------------
2844 subroutine AliHbtp_interp(rmin,rmax,rstep,nrsteps,function,
2848 CCC This routine interpolates the function values and puts the result
2849 CCC into 'answer'. It uses 2,3 or 4 mesh points which must be equally
2850 CCC spaced. The method uses the Lagrange interpolation formulas given
2851 CCC in Abramowitz and Stegun, ``Handbook of Mathematical Functions,''
2852 CCC (Dover Publications, New York, 1970); pages 878-879.
2854 CCC Definition of Variables in the Argument List:
2856 CCC rmin = lower limit of independent variable for input function
2857 CCC rmax = upper limit of independent variable for input function
2858 CCC rstep = step size of independent variable
2859 CCC nrsteps = (redundant) # of input steps
2860 CCC function(ndim) = Array of function values to be interpolated
2861 CCC ndim = array dimension size in calling program
2862 CCC r = coordinate value of independent variable where interpolation
2864 CCC answer = interpolated value
2866 CCC The algorithm will use the maximum number of points in the
2867 CCC interpolation, up to a maximum of 4
2869 CCC If the requested coordinate value, r, is out-of-range, then
2870 CCC 'answer' is returned with a 0.0 value.
2872 CCC Local Variable Type Declarations:
2874 integer*4 ndim, nrsteps, ik
2876 real*4 rmin,rmax,rstep,r,answer,rshift,p
2877 real*4 function(ndim),w1,w2,w3,w4
2881 if(abs(((rmax-rmin)/float(nrsteps-1))-rstep).gt.0.00001) then
2882 write(7,10) rmin,rmax,rstep,nrsteps
2883 10 Format(5x,'Interp mesh inconsistent:',3E12.5,I5,
2890 if(r .lt. rmin .or. r .gt. rmax) then
2891 write(7,11) rmin,rmax,r
2892 11 Format(5x,'Interp called with r out-of-range =',3E12.5)
2897 CCC Begin interpolation:
2899 if(nrsteps .eq. 2) then
2900 p = (r - rmin)/rstep
2901 answer = (1.0 - p)*function(1) + p*function(2)
2902 else if(nrsteps .eq. 3) then
2903 p = (r - (rmin + rstep))/rstep
2904 answer = 0.5*p*(p-1.0)*function(1) + (1.0 - p*p)
2905 1 *function(2) + 0.5*p*(p+1.0)*function(3)
2906 else if(nrsteps .ge. 4) then
2909 if(rshift .le. rstep) then
2911 p = (rshift - rstep)/rstep
2912 else if(rshift .ge. (rmax - rstep - rmin)) then
2914 p = (rshift - (rmax - rmin - 2.0*rstep))/rstep
2916 ik = int(rshift/rstep + 1.000001)
2917 if(ik .le. 1) ik = 2
2918 if(ik .ge. (nrsteps-1)) ik = nrsteps - 2
2919 p = (rshift - float(ik-1)*rstep)/rstep
2922 w1 = -p*(p-1.0)*(p-2.0)/6.0
2923 w2 = (p*p-1.0)*(p-2.0)/2.0
2924 w3 = -p*(p+1.0)*(p-2.0)/2.0
2925 w4 = p*(p*p-1.0)/6.0
2927 answer = w1*function(ik-1) + w2*function(ik)
2928 1 + w3*function(ik+1) + w4*function(ik+2)
2929 end if ! End # interplation points option
2934 C--------------------------------------------------------------------
2937 subroutine Hbtp_particle_prop
2940 CCC Fill particle properties table /particle/ with Geant 3 particle ID
2941 CCC numbers, charge (in units of |e|), mass in GeV/c**2 and lifetime
2942 CCC in seconds. See the Geant 3.15 Manual User's Guide, pages: CONS
2945 Include 'common_particle.inc'
2947 CCC Local Variable Type Declarations:
2951 do i = 1,part_maxlen
2955 CCC Set Particle Masses:
2957 part_mass( 1) = 0.0 ! Gamma
2958 part_mass( 2) = 0.00051099906 ! Positron
2959 part_mass( 3) = 0.00051099906 ! Electron
2960 part_mass( 4) = 0.0 ! Neutrino
2961 part_mass( 5) = 0.105658389 ! Muon+
2962 part_mass( 6) = 0.105658389 ! Muon-
2963 part_mass( 7) = 0.1349743 ! Pion0
2964 part_mass( 8) = 0.1395679 ! Pion+
2965 part_mass( 9) = 0.1395679 ! Pion-
2966 part_mass(10) = 0.497671 ! Kaon 0 long
2967 part_mass(11) = 0.493646 ! Kaon+
2968 part_mass(12) = 0.493646 ! Kaon-
2969 part_mass(13) = 0.93956563 ! Neutron
2970 part_mass(14) = 0.93827231 ! Proton
2971 part_mass(15) = 0.93827231 ! Antiproton
2972 part_mass(16) = 0.497671 ! Kaon 0 short
2973 part_mass(17) = 0.54745 ! Eta
2974 part_mass(18) = 1.11563 ! Lambda
2975 part_mass(19) = 1.18937 ! Sigma+
2976 part_mass(20) = 1.19255 ! Sigma0
2977 part_mass(21) = 1.197465 ! Sigma-
2978 part_mass(22) = 1.31485 ! Xi 0
2979 part_mass(23) = 1.32133 ! Xi -
2980 part_mass(24) = 1.67243 ! Omega
2981 part_mass(25) = 0.93956563 ! Antineutron
2982 part_mass(26) = 1.11563 ! Antilambda
2983 part_mass(27) = 1.18937 ! Anti-Sigma -
2984 part_mass(28) = 1.19255 ! Anti-Sigma 0
2985 part_mass(29) = 1.197465 ! Anti-Sigma +
2986 part_mass(30) = 1.31485 ! AntiXi 0
2987 part_mass(31) = 1.32133 ! AntiXi +
2988 part_mass(32) = 1.67243 ! Anti-Omega +
3001 part_mass(45) = 1.875613 ! Deuteron
3002 part_mass(46) = 2.80925 ! Triton
3003 part_mass(47) = 3.727417 ! Alpha
3004 part_mass(48) = 0.0 ! Geantino (Fake particle)
3005 part_mass(49) = 2.80923 ! He3
3006 part_mass(50) = 0.0 ! Cerenkov (Fake particle)
3008 CCC Set Particle Charge:
3010 part_charge( 1) = 0 ! Gamma
3011 part_charge( 2) = 1 ! Positron
3012 part_charge( 3) = -1 ! Electron
3013 part_charge( 4) = 0 ! Neutrino
3014 part_charge( 5) = 1 ! Muon+
3015 part_charge( 6) = -1 ! Muon-
3016 part_charge( 7) = 0 ! Pion0
3017 part_charge( 8) = 1 ! Pion+
3018 part_charge( 9) = -1 ! Pion-
3019 part_charge(10) = 0 ! Kaon 0 long
3020 part_charge(11) = 1 ! Kaon+
3021 part_charge(12) = -1 ! Kaon-
3022 part_charge(13) = 0 ! Neutron
3023 part_charge(14) = 1 ! Proton
3024 part_charge(15) = -1 ! Antiproton
3025 part_charge(16) = 0 ! Kaon 0 short
3026 part_charge(17) = 0 ! Eta
3027 part_charge(18) = 0 ! Lambda
3028 part_charge(19) = 1 ! Sigma+
3029 part_charge(20) = 0 ! Sigma0
3030 part_charge(21) = -1 ! Sigma-
3031 part_charge(22) = 0 ! Xi 0
3032 part_charge(23) = -1 ! Xi -
3033 part_charge(24) = -1 ! Omega
3034 part_charge(25) = 0 ! Antineutron
3035 part_charge(26) = 0 ! Antilambda
3036 part_charge(27) = -1 ! Anti-Sigma -
3037 part_charge(28) = 0 ! Anti-Sigma 0
3038 part_charge(29) = 1 ! Anti-Sigma +
3039 part_charge(30) = 0 ! AntiXi 0
3040 part_charge(31) = 1 ! AntiXi +
3041 part_charge(32) = 1 ! Anti-Omega +
3054 part_charge(45) = 1 ! Deuteron
3055 part_charge(46) = 1 ! Triton
3056 part_charge(47) = 2 ! Alpha
3057 part_charge(48) = 0 ! Geantino (Fake particle)
3058 part_charge(49) = 2 ! He3
3059 part_charge(50) = 0 ! Cerenkov (Fake particle)
3061 CCC Set Particle Lifetimes:
3063 part_lifetime( 1) = 1.0E+30 ! Gamma
3064 part_lifetime( 2) = 1.0E+30 ! Positron
3065 part_lifetime( 3) = 1.0E+30 ! Electron
3066 part_lifetime( 4) = 1.0E+30 ! Neutrino
3067 part_lifetime( 5) = 2.19703E-06 ! Muon+
3068 part_lifetime( 6) = 2.19703E-06 ! Muon-
3069 part_lifetime( 7) = 8.4E-17 ! Pion0
3070 part_lifetime( 8) = 2.603E-08 ! Pion+
3071 part_lifetime( 9) = 2.603E-08 ! Pion-
3072 part_lifetime(10) = 5.16E-08 ! Kaon 0 long
3073 part_lifetime(11) = 1.237E-08 ! Kaon+
3074 part_lifetime(12) = 1.237E-08 ! Kaon-
3075 part_lifetime(13) = 889.1 ! Neutron
3076 part_lifetime(14) = 1.0E+30 ! Proton
3077 part_lifetime(15) = 1.0E+30 ! Antiproton
3078 part_lifetime(16) = 8.922E-11 ! Kaon 0 short
3079 part_lifetime(17) = 5.53085E-19 ! Eta
3080 part_lifetime(18) = 2.632E-10 ! Lambda
3081 part_lifetime(19) = 7.99E-11 ! Sigma+
3082 part_lifetime(20) = 7.40E-20 ! Sigma0
3083 part_lifetime(21) = 1.479E-10 ! Sigma-
3084 part_lifetime(22) = 2.90E-10 ! Xi 0
3085 part_lifetime(23) = 1.639E-10 ! Xi -
3086 part_lifetime(24) = 8.22E-11 ! Omega
3087 part_lifetime(25) = 889.1 ! Antineutron
3088 part_lifetime(26) = 2.632E-10 ! Antilambda
3089 part_lifetime(27) = 7.99E-11 ! Anti-Sigma -
3090 part_lifetime(28) = 7.40E-20 ! Anti-Sigma 0
3091 part_lifetime(29) = 1.479E-10 ! Anti-Sigma +
3092 part_lifetime(30) = 2.90E-10 ! AntiXi 0
3093 part_lifetime(31) = 1.639E-10 ! AntiXi +
3094 part_lifetime(32) = 8.22E-11 ! Anti-Omega +
3095 part_lifetime(33) = 0.0
3096 part_lifetime(34) = 0.0
3097 part_lifetime(35) = 0.0
3098 part_lifetime(36) = 0.0
3099 part_lifetime(37) = 0.0
3100 part_lifetime(38) = 0.0
3101 part_lifetime(39) = 0.0
3102 part_lifetime(40) = 0.0
3103 part_lifetime(41) = 0.0
3104 part_lifetime(42) = 0.0
3105 part_lifetime(43) = 0.0
3106 part_lifetime(44) = 0.0
3107 part_lifetime(45) = 1.0E+30 ! Deuteron
3108 part_lifetime(46) = 1.0E+30 ! Triton
3109 part_lifetime(47) = 1.0E+30 ! Alpha
3110 part_lifetime(48) = 1.0E+30 ! Geantino (Fake particle)
3111 part_lifetime(49) = 1.0E+30 ! He3
3112 part_lifetime(50) = 1.0E+30 ! Cerenkov (Fake particle)
3117 C----------------------------------------------------------------
3120 subroutine correl_model
3123 CCC This subroutine computes the requested 2-body model correlation
3124 CCC function which is to be fitted by the track adjustment procedure.
3125 CCC The model values are calculated on the requested fine and coarse
3126 CCC mesh grid in momentum space. The model values are computed at the
3127 CCC mid point of each 1D bin or at the center of each 3D cell. This
3128 CCC could be refined at a later date to correspond to the integral of
3129 CCC the model function over the bin width (cell volume) divided by the
3130 CCC the bin width (cell volume).
3132 C The model includes the following options which are selected by the
3133 C 'switch*' parameters in common/parameters/:
3135 C switch_1d: 1D model as function of either Qinvar, Qtotal
3137 C switch_3d: 3D model as function of either the Bertsch-Pratt
3138 C side-out-long kinematics (but no cross term) or
3139 C the Yano-Koonin-Podgoretski perp-parallel-time
3141 C switch_type: Like and/or Unlike particles
3142 C switch_coherence: Purely incoherent source or a mixed incoherent-
3144 C switch_coulomb: Either (a) no Coulomb correction, (b) Gamow
3145 C factor, (c) NA35 parametrization, or (d) Pratt
3146 C Coulomb wave function integration for finite
3147 C size, spherical source.
3148 C switch_fermi_bose: Fermion or boson identical pairs.
3150 Include 'common_parameters.inc'
3151 Include 'common_mesh.inc'
3152 Include 'common_correlations.inc'
3154 CCC Local Variable Type Declarations:
3158 real*4 R_1dsq, Rsidesq, Routsq, Rlongsq
3159 real*4 Rperpsq, Rparallelsq, R0sq
3160 real*4 sqrtlambda,fermi_bose_sign,coulomb_factor,coherence_fac
3167 sqrtlambda = sqrt(abs(lambda))
3168 R_1dsq = ((R_1d /hbc)**2)/2.0
3169 Rsidesq = ((Rside /hbc)**2)/2.0
3170 Routsq = ((Rout /hbc)**2)/2.0
3171 Rlongsq = ((Rlong /hbc)**2)/2.0
3172 Rperpsq = ((Rperp /hbc)**2)/2.0
3173 Rparallelsq = ((Rparallel/hbc)**2)/2.0
3174 R0sq = ((R0 /hbc)**2)/2.0
3176 fermi_bose_sign = float(switch_fermi_bose)
3177 coherence_fac = switch_coherence*2.0*sqrtlambda*
3178 1 (1.0 - sqrtlambda)
3180 CCC Determine average particle pair mass for Coulomb correction:
3183 if(mass1.eq.0.0 .and. mass2.gt.0.0) massavg = mass2
3184 if(mass1.gt.0.0 .and. mass2.eq.0.0) massavg = mass1
3185 if(mass1.gt.0.0 .and. mass2.gt.0.0) massavg = 0.5*(mass1+mass2)
3187 CCC Compute 1D correlation model arrays:
3189 If(switch_1d .ge. 1) then
3190 If(n_1d_fine .gt. 0) then ! Fill the 1D fine mesh bins
3191 q = -0.5*binsize_1d_fine
3193 q = q + binsize_1d_fine
3194 b = exp(-q*q*R_1dsq)
3196 if(switch_type.eq.1 .or. switch_type.eq.3) then
3197 c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3198 1 *b*b + coherence_fac*b)
3199 if(switch_coulomb.eq.0) then
3200 coulomb_factor = 1.0
3201 else if(switch_coulomb.gt.0) then
3202 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3205 c2mod_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3208 if(switch_type.eq.2 .or. switch_type.eq.3) then
3209 c2mod_unlike_1d(i) = 1.0
3210 if(switch_coulomb.eq.0) then
3211 coulomb_factor = 1.0
3212 else if(switch_coulomb.gt.0) then
3213 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3216 c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3219 end do ! End of 1D fine mesh filling do-loop
3220 end if ! End of 1D fine mesh option
3222 If(n_1d_coarse .gt. 0) then ! Fill the 1D coarse mesh bins
3223 q = qmid_1d -0.5*binsize_1d_coarse
3224 do i = n_1d_fine + 1, n_1d_total
3225 q = q + binsize_1d_coarse
3226 b = exp(-q*q*R_1dsq)
3228 if(switch_type.eq.1 .or. switch_type.eq.3) then
3229 c2mod_like_1d(i) = 1.0 + fermi_bose_sign*(lambda
3230 1 *b*b + coherence_fac*b)
3231 if(switch_coulomb.eq.0) then
3232 coulomb_factor = 1.0
3233 else if(switch_coulomb.gt.0) then
3234 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3237 c2mod_like_1d(i) = coulomb_factor*c2mod_like_1d(i)
3240 if(switch_type.eq.2 .or. switch_type.eq.3) then
3241 c2mod_unlike_1d(i) = 1.0
3242 if(switch_coulomb.eq.0) then
3243 coulomb_factor = 1.0
3244 else if(switch_coulomb.gt.0) then
3245 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3248 c2mod_unlike_1d(i) = coulomb_factor*c2mod_unlike_1d(i)
3251 end do ! End of 1D coarse mesh filling do-loop
3252 end if ! End of 1D coarse mesh option
3253 end if ! End of 1D option
3255 CCC Compute 3D correlation model arrays:
3257 If(switch_3d .ge. 1) Then
3258 If(n_3d_fine .gt. 0) then ! Fill the 3D fine mesh bins
3259 q1 = -0.5*binsize_3d_fine
3261 q1 = q1 + binsize_3d_fine
3262 if(switch_3d.eq.1) b1=exp(-q1*q1*Rsidesq)
3263 if(switch_3d.eq.2) b1=exp(-q1*q1*Rperpsq)
3265 q2 = -0.5*binsize_3d_fine
3267 q2 = q2 + binsize_3d_fine
3268 if(switch_3d.eq.1) b2=exp(-q2*q2*Routsq)
3269 if(switch_3d.eq.2) b2=exp(-q2*q2*Rparallelsq)
3271 q3 = -0.5*binsize_3d_fine
3273 q3 = q3 + binsize_3d_fine
3274 if(switch_3d.eq.1) b3=exp(-q3*q3*Rlongsq)
3275 if(switch_3d.eq.2) b3=exp(-q3*q3*R0sq)
3278 if(switch_3d.eq.1) q = sqrt(q1*q1+q2*q2+q3*q3)
3279 if(switch_3d.eq.2) q = sqrt(q1*q1+q2*q2)
3281 if(switch_type.eq.1 .or. switch_type.eq.3) then
3282 c2mod_like_3d_fine(i,j,k) = 1.0 + fermi_bose_sign*(lambda
3283 1 *b*b + coherence_fac*b)
3284 if(switch_coulomb.eq.0) then
3285 coulomb_factor = 1.0
3286 else if(switch_coulomb.gt.0) then
3287 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3290 c2mod_like_3d_fine(i,j,k) =
3291 1 coulomb_factor*c2mod_like_3d_fine(i,j,k)
3294 if(switch_type.eq.2 .or. switch_type.eq.3) then
3295 c2mod_unlike_3d_fine(i,j,k) = 1.0
3296 if(switch_coulomb.eq.0) then
3297 coulomb_factor = 1.0
3298 else if(switch_coulomb.gt.0) then
3299 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3302 c2mod_unlike_3d_fine(i,j,k) =
3303 1 coulomb_factor*c2mod_unlike_3d_fine(i,j,k)
3308 end do ! End of 3D Fine Mesh Filling do-loops
3309 end if ! End 3D fine mesh option
3311 If(n_3d_coarse .gt. 0) then ! Fill the 3D coarse mesh bins
3312 q1 = -0.5*binsize_3d_coarse
3313 do i = 1,n_3d_coarse
3314 q1 = q1 + binsize_3d_coarse
3315 if(switch_3d.eq.1) b1=exp(-q1*q1*Rsidesq)
3316 if(switch_3d.eq.2) b1=exp(-q1*q1*Rperpsq)
3318 q2 = -0.5*binsize_3d_coarse
3319 do j = 1,n_3d_coarse
3320 q2 = q2 + binsize_3d_coarse
3321 if(switch_3d.eq.1) b2=exp(-q2*q2*Routsq)
3322 if(switch_3d.eq.2) b2=exp(-q2*q2*Rparallelsq)
3324 q3 = -0.5*binsize_3d_coarse
3325 do k = 1,n_3d_coarse
3326 q3 = q3 + binsize_3d_coarse
3327 if(switch_3d.eq.1) b3=exp(-q3*q3*Rlongsq)
3328 if(switch_3d.eq.2) b3=exp(-q3*q3*R0sq)
3331 if(switch_3d.eq.1) q = sqrt(q1*q1+q2*q2+q3*q3)
3332 if(switch_3d.eq.2) q = sqrt(q1*q1+q2*q2)
3334 if(switch_type.eq.1 .or. switch_type.eq.3) then
3335 c2mod_like_3d_coarse(i,j,k) = 1.0+fermi_bose_sign*(lambda
3336 1 *b*b + coherence_fac*b)
3337 if(switch_coulomb.eq.0) then
3338 coulomb_factor = 1.0
3339 else if(switch_coulomb.gt.0) then
3340 Call coulomb(switch_coulomb,q,1,massavg,Q0,
3343 c2mod_like_3d_coarse(i,j,k) =
3344 1 coulomb_factor*c2mod_like_3d_coarse(i,j,k)
3347 if(switch_type.eq.2 .or. switch_type.eq.3) then
3348 c2mod_unlike_3d_coarse(i,j,k) = 1.0
3349 if(switch_coulomb.eq.0) then
3350 coulomb_factor = 1.0
3351 else if(switch_coulomb.gt.0) then
3352 Call coulomb(switch_coulomb,q,-1,massavg,Q0,
3355 c2mod_unlike_3d_coarse(i,j,k) =
3356 1 coulomb_factor*c2mod_unlike_3d_coarse(i,j,k)
3361 end do ! End of 3D Coarse Mesh Filling do-loops
3362 c2mod_like_3d_coarse(1,1,1) = 0.0
3363 c2mod_unlike_3d_coarse(1,1,1) = 0.0
3364 end if ! End of 3D Coarse Mesh Option
3366 End If ! End of 3D Option
3371 C----------------------------------------------------------------------
3374 subroutine coulomb(control,q,sign,mass,Q0,factor)
3377 CCC Compute Coulomb correction to the two-body correlation functions
3378 C for like and unlike charges for particles of the same mass.
3379 C Three methods are allowed:
3381 C If control = 1, Then use the gamow factor for point sources
3382 C If control = 2, Then use the NA35 finite source size empirical
3383 C correction factor from eq.(5) in Z. Phys. C73,
3385 C If control = 3, Then use the Pratt finite source size, numerically
3386 C integrated Coulomb correction factor with inter-
3389 C Other parameters in the argument list are:
3391 C q = 3-vector momentum difference for track pair, in GeV/c.
3392 C sign = algebraic sign of the charge product for the track pair.
3393 C mass = particle mass in GeV, it is assumed that both particles
3394 C have the same mass, e.g. pi+ and pi-, but not K+ and pi-.
3395 C Q0 = NA35 parameter in GeV/c if control = 2
3396 C = Source radius in fm if control = 3
3397 C factor = Multiplicative Coulomb correction result which is
3398 C calculated here and returned to the calling program.
3401 Include 'common_coulomb.inc'
3403 CCC Local Variable Type Declarations:
3405 integer*4 control, sign
3407 real*4 pi,q,mass,Q0,factor,alpha,eta,eta2pi
3409 parameter (pi = 3.141592654)
3410 parameter (alpha = 0.00729735)
3412 CCC Compute Gamow factor for control options 1 and 2:
3414 if(control.eq.1 .or. control.eq.2) then
3415 if(q .le. 0.001) then
3416 if(sign .gt. 0) gamow = 0.0
3417 if(sign .lt. 0) gamow = 86.0
3419 eta = sign*mass*alpha/q
3421 gamow = eta2pi/(exp(eta2pi) - 1.0)
3425 CCC Compute Coulomb Correction factor for options 1, 2 and 3:
3427 if(control .eq. 1) then
3429 else if(control .eq. 2) then
3430 factor = 1.0 + (gamow - 1.0)*exp(-q/Q0)
3431 else if(control .eq. 3) then
3433 if(q .le. q_coul(1)) then
3434 if(sign .gt. 0) factor = c2_coul_like(1)
3435 if(sign .lt. 0) factor = c2_coul_unlike(1)
3436 else if(q .ge. q_coul(max_c2_coul - 1)) then
3437 if(sign .gt. 0) factor = c2_coul_like(max_c2_coul - 1)
3438 if(sign .lt. 0) factor = c2_coul_unlike(max_c2_coul - 1)
3440 if(sign .gt. 0) then
3441 Call LAGRNG1(q,q_coul,factor,c2_coul_like,
3442 1 max_c2_coul,1,5,max_c2_coul,1)
3443 else if(sign .lt. 0) then
3444 Call LAGRNG1(q,q_coul,factor,c2_coul_unlike,
3445 1 max_c2_coul,1,5,max_c2_coul,1)
3448 end if ! END Coulomb correction evaluation, control selection opt.
3453 C---------------------------------------------------------------------
3456 SUBROUTINE LAGRNG1 (X,ARG,Y,VAL,NDIM,NFS,NPTS,MAXARG,MAXFS)
3457 IMPLICIT REAL*4(A-H,O-Z)
3459 C LAGRANGE INTERPOLATION,UNEQUALLY SPACED POINTS
3460 C ROUTINE OBTAINED FROM R. LANDAU, UNIV. OF OREGON.
3461 C ARG=VECTOR OF INDEPENDENT VARIABLE CONTAINING MAXARG VALUES.
3462 C VAL=MATRIX OF FUNCTION VALUES CORRESPONDING TO ARG. (MAXFS
3463 C FUNCTIONS AT MAXARG VALUES.)
3464 C X =VALUE OF INDEP. VARIABLE FOR WHICH INTERPOLATION IS DESIRED.
3465 C Y =VECTOR OF MAXFS FUNCTION VALUES RESULTING FROM SIMUL. INTERP.
3466 C NDIM=NUMBER OF ARG VALUES TO BE USED. (NDIM.LE.MAXARG)
3467 C NFS=NUMBER OF FUNCTIONS SIMUL. INTERP (NFS.LE.MAXFS)
3468 C NPTS=NUMBER OF POINTS USED IN INTERPOLATION. (NPTS=2,3,4,5,6)
3470 DIMENSION ARG(MAXARG), VAL(MAXFS,MAXARG), Y(MAXFS)
3472 C -----FIND X0, THE CLOSEST POINT TO X.
3476 10 IF ((X.LE.ARG(NI)).OR.(X.GE.ARG(NF))) GO TO 30
3477 IF ((NF-NI+1).EQ.2) GO TO 70
3479 IF (X.GT.ARG(NMID)) GO TO 20
3485 C ------ X IS ONE OF THE TABLULATED VALUES.
3487 30 IF (X.LE.ARG(NI)) GO TO 60
3496 C ------- 2 PTS LEFT, CHOOSE SMALLER ONE.
3500 GO TO (110,100,90,80), NN
3502 IF (((N0+3).GT.NDIM).OR.((N0-2).LT.1)) GO TO 90
3506 IF ((N0+2).GT.NDIM) GO TO 110
3507 IF ((N0-2).LT.1) GO TO 100
3511 IF (((N0+2).GT.NDIM).OR.((N0-1).LT.1)) GO TO 110
3514 110 IF ((N0+1).LT.NDIM) GO TO 120
3516 C ------N0=NDIM, SPECIAL CASE.
3521 IF ((N0-1).LT.1) NUSED=2
3524 C ------AT LEAST 2 PTS LEFT.
3531 IF (NUSED.EQ.2) GO TO 140
3533 C ------AT LEAST 3 PTS.
3538 CM1=-Y0*Y1/Y0M1/YM11
3541 IF (NUSED.EQ.3) GO TO 160
3543 C ------AT LEAST 4 PTS
3552 C2=-YM1*Y0*Y1/YM12/Y02/Y12
3553 IF (NUSED.EQ.4) GO TO 180
3555 C ------AT LEAST 5 PTS.
3562 CM2=YM1*Y0*Y1*Y2/YM2M1/YM20/YM21/YM22
3567 IF (NUSED.EQ.5) GO TO 200
3569 C ------AT LEAST 6 PTS.
3582 C3=YM2*YM1*Y0*Y1*Y2/YM23/YM13/Y03/Y13/Y23
3586 150 Y(N)=C0*VAL(N,N0)+C1*VAL(N,N0+1)
3590 170 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)
3594 190 Y(N)=CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C2*VAL(N,N0+2)
3598 210 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
3603 230 Y(N)=CM2*VAL(N,N0-2)+CM1*VAL(N,N0-1)+C0*VAL(N,N0)+C1*VAL(N,N0+1)+C
3604 12*VAL(N,N0+2)+C3*VAL(N,N0+3)
3609 C-------------------------------------------------------------------
3612 subroutine Hbtp_kin(px,py,pz,E,pt,phi,eta,mass,control)
3615 CCC Four-momentum kinematics conversion:
3617 C If control = 1, use input {px,py,pz,mass} to calculate
3619 C If control = 2, use input {pt,phi,eta,mass} to calculate
3622 C Units: Momentum are in GeV/c
3623 C Energy and mass are in GeV
3624 C Angles are in degrees
3626 CCC Local Variable Type Declarations:
3630 real*4 px,py,pz,E,pt,phi,eta,mass
3631 real*4 theta,pi,rad,pcut,x,y
3632 parameter (pi = 3.141592654)
3633 parameter (pcut = 0.000001)
3637 If(control .eq. 1) Then ! Use {px,py,pz,mass} --> {E,pt,phi,eta}
3638 pt = sqrt(px*px + py*py)
3639 E = sqrt(pt*pt + pz*pz + mass*mass)
3641 CCC Compute azimuthal angle phi; treat pt = 0.0 and py = 0.0 cases
3644 if(pt .le. pcut) then
3646 else if(pt.gt.pcut.and.abs(py).le.pcut.and.px.lt.0.0) then
3651 if(phi .lt. 0.0) phi = phi + 2.0*pi
3654 CCC Compute pseudorapidity:
3656 if(pt.le.pcut .and. abs(pz).le.pcut) then
3658 else if(pt.le.pcut .and. abs(pz).gt.pcut) then
3659 eta = 0.5*log((E+pz)/(E-pz)) ! Use beam rapidity
3661 theta = atan2(pt,pz)
3662 eta = -log(tan(theta/2.0))
3665 Else If(control .eq. 2) Then ! Use {pt,phi,eta,mass} --> {E,px,py,pz}
3667 px = pt*cos(phi/rad)
3668 py = pt*sin(phi/rad)
3669 if(abs(eta) .le. pcut) then
3675 C theta = 2.0*atan(exp(-eta))
3679 E = sqrt(pt*pt + pz*pz + mass*mass)
3681 End If ! End control options
3686 C----------------------------------------------------------------------
3689 subroutine qdiff(px1,py1,pz1,E1,px2,py2,pz2,E2,qinvar2,qtotal2,
3690 1 qvector2,qside2,qout2,qlong2,qperp2,qtime2)
3693 CCC This subroutine computes the various relative momenta for given
3694 CCC input 4-momentum for particles 1 and 2. The subr: returns the
3695 CCC square of the momentum. All energy and momenta are in GeV.
3697 C Input 4-momentum for particle 1: {px1,py1,pz1,E1}
3698 C Input 4-momentum for particle 2: {px2,py2,pz2,E2}
3700 CCC Computed Momentum Differences are the following:
3702 C qinvar2 = Q-invariant**2
3703 C qtotal2 = Q-Total**2 (space**2 + time**2)
3704 C qvector2 = 3-momentum vector difference squared
3705 C qside2 = q-side**2 of Bertsch-Pratt 3D source models
3706 C qout2 = q-out**2 of Bertsch-Pratt 3D source models
3707 C qlong2 = q-long**2 of Bertsch-Pratt 3D source models
3708 C qperp2 = q-perpendicular**2 of YKP 3D source models
3709 C qparallel2= q-parallel**2 of YKP 3D source models
3710 C = qlong2 (not assigned a separate variable name)
3711 C qtime2 = q-time-like**2 of YKP 3D source models
3713 CCC Local Variable Type Declarations:
3715 real*4 px1,py1,pz1,E1,px2,py2,pz2,E2
3716 real*4 qinvar2,qtotal2,qvector2
3717 real*4 qside2,qout2,qlong2
3718 real*4 qperp2,qtime2
3719 real*4 px12sq,py12sq,pz12sq,E12sq
3721 px12sq = (px1 - px2)**2
3722 py12sq = (py1 - py2)**2
3723 pz12sq = (pz1 - pz2)**2
3724 E12sq = (E1 - E2)**2
3726 qvector2 = px12sq + py12sq + pz12sq
3727 qinvar2 = qvector2 - E12sq
3728 qtotal2 = qvector2 + E12sq
3731 qout2 = (px1*px1 - px2*px2 + py1*py1 - py2*py2)
3732 qout2 = qout2*qout2/((px1+px2)**2 + (py1+py2)**2)
3733 qperp2 = px12sq + py12sq
3734 qside2 = qperp2 - qout2
3736 if (qside2 .lt. 0) then
3737 C write(*,*) 'qside2 is less then 0', qside2
3738 C write(*,*) ' qperp2, qout2', qperp2, qout2
3739 C write(*,*) ' px1,py1,pz1,E1 ',px1,py1,pz1,E1
3740 C write(*,*) ' px2,py2,pz2,E2 ',px2,py2,pz2,E2
3746 C-----------------------------------------------------------------------
3749 subroutine mean_rms(a,ndim,npts,mean,rms)
3752 CCC Calculate the mean and standard deviation (rms) for input
3753 C distribution a() for npts number of values.
3754 C ndim = dimension of array a() in calling program.
3756 CCC Local Variable Type Declarations:
3758 integer*4 ndim, npts, i
3759 real*4 a(ndim), mean, rms, sum_mean, sum_rms
3761 if(npts .le. 0) then
3765 else if(npts .eq. 1) then
3773 sum_mean = sum_mean + a(i)
3775 mean = sum_mean/float(npts)
3778 sum_rms = sum_rms + (a(i) - mean)**2
3780 rms = sqrt((sum_rms)/float(npts - 1))
3786 C-----------------------------------------------------------------------
3789 subroutine tindex(mode,track_id)
3792 CCC This subroutine locates tracks in {px,py,pz} sectors
3793 C and sets the sector index numbers
3794 C in track table 'trk' and/or 'trk2', depending on the value of 'mode'.
3795 C If a track's momentum is out of the sector ranges, then the track
3796 C will be assigned to, and counted in the nearest sector cell on the
3799 C If mode = 1, apply to tracks in table 'trk'
3800 C If mode = 2, apply to tracks in table 'trk2'
3802 C If track_id = 0, then do this for all tracks
3803 C If track_id = i (where i.gt.0) then do this for track row i only
3805 Include 'common_parameters.inc'
3806 Include 'common_mesh.inc'
3808 Include 'common_track.inc'
3809 Include 'common_track2.inc'
3811 CCC Local Variable Type Declarations:
3813 integer*4 i,track_id,mode
3815 C-----------------------
3817 C-----------------------
3819 If(track_id .eq. 0) Then
3820 do i = 1,n_part_tot_trk
3821 trk_px_sec(i) = int(((trk_px(i) - px_min)/delpx)+1.00001)
3822 if(trk_px_sec(i) .lt.1) trk_px_sec(i) = 1
3823 if(trk_px_sec(i) .gt. n_px_bins) trk_px_sec(i) = n_px_bins
3824 trk_py_sec(i) = int(((trk_py(i) - py_min)/delpy)+1.00001)
3825 if(trk_py_sec(i) .lt.1) trk_py_sec(i) = 1
3826 if(trk_py_sec(i) .gt. n_py_bins) trk_py_sec(i) = n_py_bins
3827 trk_pz_sec(i) = int(((trk_pz(i) - pz_min)/delpz)+1.00001)
3828 if(trk_pz_sec(i) .lt.1) trk_pz_sec(i) = 1
3829 if(trk_pz_sec(i) .gt. n_pz_bins) trk_pz_sec(i) = n_pz_bins
3830 trk_sector(i) = trk_px_sec(i) + (trk_py_sec(i) - 1)*
3831 1 n_px_bins + (trk_pz_sec(i) - 1)*n_px_bins*n_py_bins
3834 Else If(track_id .gt. 0) Then
3836 trk_px_sec(i) = int(((trk_px(i) - px_min)/delpx)+1.00001)
3837 if(trk_px_sec(i) .lt.1) trk_px_sec(i) = 1
3838 if(trk_px_sec(i) .gt. n_px_bins) trk_px_sec(i) = n_px_bins
3839 trk_py_sec(i) = int(((trk_py(i) - py_min)/delpy)+1.00001)
3840 if(trk_py_sec(i) .lt.1) trk_py_sec(i) = 1
3841 if(trk_py_sec(i) .gt. n_py_bins) trk_py_sec(i) = n_py_bins
3842 trk_pz_sec(i) = int(((trk_pz(i) - pz_min)/delpz)+1.00001)
3843 if(trk_pz_sec(i) .lt.1) trk_pz_sec(i) = 1
3844 if(trk_pz_sec(i) .gt. n_pz_bins) trk_pz_sec(i) = n_pz_bins
3845 trk_sector(i) = trk_px_sec(i) + (trk_py_sec(i) - 1)*
3846 1 n_px_bins + (trk_pz_sec(i) - 1)*n_px_bins*n_py_bins
3849 C-----------------------------
3850 Else If (mode.eq.2) Then
3851 C-----------------------------
3853 If(track_id .eq. 0) Then
3854 do i = 1,n_part_tot_trk2
3855 trk2_px_sec(i) = int(((trk2_px(i) - px_min)/delpx)+1.00001)
3856 if(trk2_px_sec(i) .lt.1) trk2_px_sec(i) = 1
3857 if(trk2_px_sec(i) .gt. n_px_bins) trk2_px_sec(i) = n_px_bins
3858 trk2_py_sec(i) = int(((trk2_py(i) - py_min)/delpy)+1.00001)
3859 if(trk2_py_sec(i) .lt.1) trk2_py_sec(i) = 1
3860 if(trk2_py_sec(i) .gt. n_py_bins) trk2_py_sec(i) = n_py_bins
3861 trk2_pz_sec(i) = int(((trk2_pz(i) - pz_min)/delpz)+1.00001)
3862 if(trk2_pz_sec(i) .lt.1) trk2_pz_sec(i) = 1
3863 if(trk2_pz_sec(i) .gt. n_pz_bins) trk2_pz_sec(i) = n_pz_bins
3864 trk2_sector(i) = trk2_px_sec(i) + (trk2_py_sec(i) - 1)*
3865 1 n_px_bins + (trk2_pz_sec(i) - 1)*n_px_bins*n_py_bins
3868 Else If(track_id .gt. 0) Then
3870 trk2_px_sec(i) = int(((trk2_px(i) - px_min)/delpx)+1.00001)
3871 if(trk2_px_sec(i) .lt.1) trk2_px_sec(i) = 1
3872 if(trk2_px_sec(i) .gt. n_px_bins) trk2_px_sec(i) = n_px_bins
3873 trk2_py_sec(i) = int(((trk2_py(i) - py_min)/delpy)+1.00001)
3874 if(trk2_py_sec(i) .lt.1) trk2_py_sec(i) = 1
3875 if(trk2_py_sec(i) .gt. n_py_bins) trk2_py_sec(i) = n_py_bins
3876 trk2_pz_sec(i) = int(((trk2_pz(i) - pz_min)/delpz)+1.00001)
3877 if(trk2_pz_sec(i) .lt.1) trk2_pz_sec(i) = 1
3878 if(trk2_pz_sec(i) .gt. n_pz_bins) trk2_pz_sec(i) = n_pz_bins
3879 trk2_sector(i) = trk2_px_sec(i) + (trk2_py_sec(i) - 1)*
3880 1 n_px_bins + (trk2_pz_sec(i) - 1)*n_px_bins*n_py_bins
3890 C------------------------------------------------------------------------
3893 subroutine stm_build(mode,track_index,old_sector)
3896 CCC This subroutine fills or updates the track-sector information
3897 C table sec_trk_map or, for the reference calculations, it fills
3898 C sec_trk_map2. These track-sector tables contain the information
3899 C about the track occupancy, status, etc. for all of the {px,py,pz}
3904 C If track_index = 0, then fill information for all tracks in 'trk',
3906 C If track_index = i (where i.gt.0) then fill only the track-sector
3907 C information for track i in 'trk', into table 'stm'
3908 C Also for this case, if old_sector .ne. 0, then
3909 C remove the track # and ID information for this
3910 C old sector # from table stm
3912 C For Mode = 2: Fill information for all tracks in 'trk2' into table 'stm2'
3914 Include 'common_parameters.inc'
3915 Include 'common_mesh.inc'
3917 Include 'common_track.inc'
3918 Include 'common_track2.inc'
3919 Include 'common_sec_track.inc'
3920 Include 'common_sec_track2.inc'
3922 CCC Local Variable Type Declarations:
3924 integer*4 i,j,mode,track_index,old_sector,row
3925 integer*4 temp(max_trk_sec)
3927 C------------------------
3929 C------------------------
3931 If (track_index .eq. 0) Then
3934 stm_n_trk_sec(i) = 0
3935 do j = 1,max_trk_sec
3936 stm_track_id(j,i) = 0
3943 do i = 1,n_part_tot_trk
3944 if(trk_flag(i) .eq. 0) then
3946 stm_n_trk_sec(row) = stm_n_trk_sec(row) + 1
3947 if(stm_n_trk_sec(row) .le. max_trk_sec) then
3948 stm_track_id(stm_n_trk_sec(row),row) = trk_id(i)
3952 stm_n_trk_sec(row) = stm_n_trk_sec(row) - 1
3960 Else If (track_index .gt. 0) Then
3962 if(old_sector .ne. 0) then
3963 CCC Remove track from old sector:
3965 do i = 1,stm_n_trk_sec(old_sector)
3966 if(stm_track_id(i,old_sector) .ne. track_index) then
3968 temp(j) = stm_track_id(i,old_sector)
3971 stm_n_trk_sec(old_sector) = j
3972 do i = 1,max_trk_sec
3973 stm_track_id(i,old_sector) = 0
3975 do i = 1,stm_n_trk_sec(old_sector)
3976 stm_track_id(i,old_sector) = temp(i)
3979 CCC Update with new sector location of track:
3981 if(trk_flag(i) .eq. 0) then
3983 stm_n_trk_sec(row) = stm_n_trk_sec(row) + 1
3984 if(stm_n_trk_sec(row) .le. max_trk_sec) then
3985 stm_track_id(stm_n_trk_sec(row),row) = trk_id(i)
3989 stm_n_trk_sec(row) = stm_n_trk_sec(row) - 1
3997 C-----------------------------
3998 Else If (mode.eq.2) Then
3999 C-----------------------------
4001 If (track_index .eq. 0) Then
4002 do i = 1,sec_maxlen2
4004 stm2_n_trk_sec(i) = 0
4005 do j = 1,max_trk_sec2
4006 stm2_track_id(j,i) = 0
4013 do i = 1,n_part_tot_trk2
4014 if(trk2_flag(i) .eq. 0) then
4015 row = trk2_sector(i)
4016 stm2_n_trk_sec(row) = stm2_n_trk_sec(row) + 1
4017 if(stm2_n_trk_sec(row) .le. max_trk_sec2) then
4018 stm2_track_id(stm2_n_trk_sec(row),row) = trk2_id(i)
4022 stm2_n_trk_sec(row) = stm2_n_trk_sec(row) - 1
4032 End If ! End mode = 1,2 selection options
4038 C-----------------------------------------------------------------------
4041 subroutine sec_index(index,nbins,index_min,index_max)
4044 CCC Calculate track-sector neighboring bins and min->max range:
4046 CCC Local Variable Type Declarations:
4048 integer*4 index,nbins,index_min,index_max
4050 index_min = index - 1
4051 if(index_min .lt. 1) index_min = 1
4052 index_max = index + 1
4053 if(index_max .gt. nbins) index_max = nbins
4058 C-----------------------------------------------------------------------
4061 subroutine dist_range(mode,ntracks_out,ntracks_flagged)
4064 CCC Determine if tracks are out of acceptance range in pt, phi and eta,
4065 C and, if so, then set the 'out_flag' variable in the track table 'trk'
4067 CCC For Mode = 1, use track table 'trk'
4068 CCC For Mode = 2, use track table 'trk2'
4070 CCC Count the number of flagged tracks, i.e. trk(i).flag = 1, and "out
4071 C of acceptance range" tracks, i.e. trk(i).out_flag = 1, for both
4072 C particle ID types. Determine the number of tracks to use in the
4073 C correlation fit for each particle ID type.
4075 Include 'common_parameters.inc'
4076 Include 'common_mesh.inc'
4078 Include 'common_track.inc'
4079 Include 'common_track2.inc'
4081 CCC Local Variable Type Declarations:
4083 integer*4 i,mode,ntracks_out,ntracks_flagged
4085 C------------------------
4087 C------------------------
4095 do i = 1,n_part_tot_trk
4096 if(trk_flag(i) .eq. 0) then
4097 if(trk_pt(i) .lt. pt_min .or. trk_pt(i) .gt. pt_max)
4099 if(trk_phi(i).lt.phi_min .or. trk_phi(i).gt.phi_max)
4101 if(trk_eta(i).lt.eta_min .or. trk_eta(i).gt.eta_max)
4103 else if(trk_flag(i) .eq. 1) then
4104 ntracks_flagged = ntracks_flagged + 1
4109 do i = 1,n_part_tot_trk
4110 if(trk_out_flag(i) .eq. 1) ntracks_out = ntracks_out + 1
4113 n_part_used_1_trk = 0
4114 n_part_used_2_trk = 0
4115 do i = 1,n_part_tot_trk
4116 if(trk_flag(i) .eq. 0) then
4117 if(trk_ge_pid(i) .eq. pid(1)) then
4118 n_part_used_1_trk = n_part_used_1_trk + 1
4119 else if(trk_ge_pid(i) .eq. pid(2)) then
4120 n_part_used_2_trk = n_part_used_2_trk + 1
4125 C-----------------------------
4126 Else If (mode.eq.2) Then
4127 C-----------------------------
4129 do i = 1,trk2_maxlen
4130 trk2_out_flag(i) = 0
4135 do i = 1,n_part_tot_trk2
4136 if(trk2_flag(i) .eq. 0) then
4137 if(trk2_pt(i) .lt. pt_min .or. trk2_pt(i) .gt. pt_max)
4138 1 trk2_out_flag(i)=1
4139 if(trk2_phi(i).lt.phi_min .or. trk2_phi(i).gt.phi_max)
4140 1 trk2_out_flag(i)=1
4141 if(trk2_eta(i).lt.eta_min .or. trk2_eta(i).gt.eta_max)
4142 1 trk2_out_flag(i)=1
4143 else if(trk2_flag(i) .eq. 1) then
4144 ntracks_flagged = ntracks_flagged + 1
4149 do i = 1,n_part_tot_trk2
4150 if(trk2_out_flag(i) .eq. 1) ntracks_out = ntracks_out + 1
4153 n_part_used_1_trk2 = 0
4154 n_part_used_2_trk2 = 0
4155 do i = 1,n_part_tot_trk2
4156 if(trk2_flag(i) .eq. 0) then
4157 if(trk2_ge_pid(i) .eq. pid(1)) then
4158 n_part_used_1_trk2 = n_part_used_1_trk2 + 1
4159 else if(trk2_ge_pid(i) .eq. pid(2)) then
4160 n_part_used_2_trk2 = n_part_used_2_trk2 + 1
4166 End If ! End mode = 1,2 selection option
4172 C--------------------------------------------------------------------
4175 subroutine histog1(mode,itrack,pid_index,pidnum,pt_save,
4176 1 phi_save,eta_save)
4179 CCC This subroutine computes and/or updates the one-body histograms
4180 C according to the selected 'mode' and for the requested particle
4183 C Note: If the track momentum is out-of-range in {pt,phi,eta},
4184 C then it is ignored. The {pt,phi,eta} dependences for
4185 C the 1-dimensional histogramming are treated independently.
4186 C It is therefore possible for the sum of the number of
4187 C particles in the pt, phi and eta one-body, 1D histograms
4190 CCC Mode = 1, Fill the one-body histograms (hist1*) for selected
4191 C particle ID type, for the initial input distribution,
4192 C using the momenta in 'trk'
4194 C Mode = 2, Remove particle 'itrack' from temporary one-body hist-
4195 C ogram (htmp1*) for selected particle ID type, using
4196 C momentum values given by pt_save, phi_save, eta_save.
4198 C Mode = 3, Add particle 'itrack' to the temporary one-body hist-
4199 C ogram (htmp1*) for selected particle ID type, using
4200 C momentum values in track table 'trk'.
4202 C Mode = 4, Fill the one-body histograms (hist1*) for selected
4203 C particle ID type, for the initial input distribution,
4204 C using the momenta in 'trk2'
4206 C itrack = track index # for the removed or added track for mode =
4207 C 2 or 3, respectively.
4209 C pid_index = 1 or 2 for the first or second particle ID type, and
4210 C for filling/update either hist1*1 or hist1*2, and
4211 C similarly for htmp1*1 or htmp1*2
4213 C pidnum = Geant particle ID # for the track(s) to be filled or
4216 C pt_save = Removed track's pt value.
4218 C phi_save = Removed track's phi value.
4220 C eta_save = Removed track's eta value.
4222 Include 'common_parameters.inc'
4223 Include 'common_mesh.inc'
4224 Include 'common_histograms.inc'
4226 Include 'common_track.inc'
4227 Include 'common_track2.inc'
4229 CCC Local Variable Type Declarations:
4231 integer*4 mode,itrack,i,pid_index,pidnum,index
4232 integer*4 trk_counter,trk2_counter
4234 real*4 pt_save,phi_save,eta_save
4235 real*4 delpt,delphi,deleta
4237 C-------------------------
4239 C-------------------------
4241 CCC Fill one-body histograms for requested particle ID from table 'trk'
4242 CCC Initialize necessary arrays to zero:
4244 if(pid_index .eq. 1) then
4250 else if(pid_index .eq. 2) then
4260 do i = 1,n_part_tot_trk
4261 if(trk_ge_pid(i).eq.pidnum .and. trk_flag(i).eq.0) then
4262 trk_counter = trk_counter + 1
4263 delpt = trk_pt(i) - pt_min
4264 delphi = trk_phi(i) - phi_min
4265 deleta = trk_eta(i) - eta_min
4267 index = int((delpt/pt_bin_size) + 0.99999)
4268 if(index.ge.1 .and. index.le.max_h_1d) then
4269 if(pid_index.eq.1) hist1_pt_1(index) =
4270 1 hist1_pt_1(index) + 1
4271 if(pid_index.eq.2) hist1_pt_2(index) =
4272 1 hist1_pt_2(index) + 1
4275 index = int((delphi/phi_bin_size) + 0.99999)
4276 if(index.ge.1 .and. index.le.max_h_1d) then
4277 if(pid_index.eq.1) hist1_phi_1(index) =
4278 1 hist1_phi_1(index) + 1
4279 if(pid_index.eq.2) hist1_phi_2(index) =
4280 1 hist1_phi_2(index) + 1
4283 index = int((deleta/eta_bin_size) + 0.99999)
4284 if(index.ge.1 .and. index.le.max_h_1d) then
4285 if(pid_index.eq.1) hist1_eta_1(index) =
4286 1 hist1_eta_1(index) + 1
4287 if(pid_index.eq.2) hist1_eta_2(index) =
4288 1 hist1_eta_2(index) + 1
4294 if(pid_index .eq. 1) n_part_used_1_trk = trk_counter
4295 if(pid_index .eq. 2) n_part_used_2_trk = trk_counter
4297 C--------------------------------
4298 Else If (mode .eq. 2) Then
4299 C--------------------------------
4301 CCC Remove track # 'itrack' from fitting histograms in htmp1*,
4302 CCC use pt_save, phi_save, eta_save for the old momentum values
4303 CCC in order to determine which bins to decrement.
4305 if(trk_ge_pid(itrack).eq.pidnum.and.trk_flag(itrack).eq.0)then
4306 delpt = pt_save - pt_min
4307 delphi = phi_save - phi_min
4308 deleta = eta_save - eta_min
4310 index = int((delpt/pt_bin_size) + 0.99999)
4311 if(index.ge.1 .and. index.le.max_h_1d) then
4312 if(pid_index.eq.1) htmp1_pt_1(index) =
4313 1 htmp1_pt_1(index) - 1
4314 if(pid_index.eq.2) htmp1_pt_2(index) =
4315 1 htmp1_pt_2(index) - 1
4318 index = int((delphi/phi_bin_size) + 0.99999)
4319 if(index.ge.1 .and. index.le.max_h_1d) then
4320 if(pid_index.eq.1) htmp1_phi_1(index) =
4321 1 htmp1_phi_1(index) - 1
4322 if(pid_index.eq.2) htmp1_phi_2(index) =
4323 1 htmp1_phi_2(index) - 1
4326 index = int((deleta/eta_bin_size) + 0.99999)
4327 if(index.ge.1 .and. index.le.max_h_1d) then
4328 if(pid_index.eq.1) htmp1_eta_1(index) =
4329 1 htmp1_eta_1(index) - 1
4330 if(pid_index.eq.2) htmp1_eta_2(index) =
4331 1 htmp1_eta_2(index) - 1
4336 C--------------------------------
4337 Else If (mode .eq. 3) Then
4338 C--------------------------------
4340 CCC Add track # 'itrack' to fitting histograms in htmp1*,
4341 CCC use pt, phi and eta values in track table 'trk(itrack)'
4342 CCC for the new/added track position.
4344 if(trk_ge_pid(itrack).eq.pidnum.and.trk_flag(itrack).eq.0)then
4345 delpt = trk_pt(itrack) - pt_min
4346 delphi = trk_phi(itrack) - phi_min
4347 deleta = trk_eta(itrack) - eta_min
4349 index = int((delpt/pt_bin_size) + 0.99999)
4350 if(index.ge.1 .and. index.le.max_h_1d) then
4351 if(pid_index.eq.1) htmp1_pt_1(index) =
4352 1 htmp1_pt_1(index) + 1
4353 if(pid_index.eq.2) htmp1_pt_2(index) =
4354 1 htmp1_pt_2(index) + 1
4357 index = int((delphi/phi_bin_size) + 0.99999)
4358 if(index.ge.1 .and. index.le.max_h_1d) then
4359 if(pid_index.eq.1) htmp1_phi_1(index) =
4360 1 htmp1_phi_1(index) + 1
4361 if(pid_index.eq.2) htmp1_phi_2(index) =
4362 1 htmp1_phi_2(index) + 1
4365 index = int((deleta/eta_bin_size) + 0.99999)
4366 if(index.ge.1 .and. index.le.max_h_1d) then
4367 if(pid_index.eq.1) htmp1_eta_1(index) =
4368 1 htmp1_eta_1(index) + 1
4369 if(pid_index.eq.2) htmp1_eta_2(index) =
4370 1 htmp1_eta_2(index) + 1
4375 C------------------------------
4376 Else If (mode.eq.4) Then
4377 C------------------------------
4379 CCC Fill one-body histograms for requested particle ID from table 'trk2'
4380 CCC Initialize necessary arrays to zero:
4382 if(pid_index .eq. 1) then
4388 else if(pid_index .eq. 2) then
4398 do i = 1,n_part_tot_trk2
4399 if(trk2_ge_pid(i).eq.pidnum .and. trk2_flag(i).eq.0) then
4400 trk2_counter = trk2_counter + 1
4401 delpt = trk2_pt(i) - pt_min
4402 delphi = trk2_phi(i) - phi_min
4403 deleta = trk2_eta(i) - eta_min
4405 index = int((delpt/pt_bin_size) + 0.99999)
4406 if(index.ge.1 .and. index.le.max_h_1d) then
4407 if(pid_index.eq.1) hist1_pt_1(index) =
4408 1 hist1_pt_1(index) + 1
4409 if(pid_index.eq.2) hist1_pt_2(index) =
4410 1 hist1_pt_2(index) + 1
4413 index = int((delphi/phi_bin_size) + 0.99999)
4414 if(index.ge.1 .and. index.le.max_h_1d) then
4415 if(pid_index.eq.1) hist1_phi_1(index) =
4416 1 hist1_phi_1(index) + 1
4417 if(pid_index.eq.2) hist1_phi_2(index) =
4418 1 hist1_phi_2(index) + 1
4421 index = int((deleta/eta_bin_size) + 0.99999)
4422 if(index.ge.1 .and. index.le.max_h_1d) then
4423 if(pid_index.eq.1) hist1_eta_1(index) =
4424 1 hist1_eta_1(index) + 1
4425 if(pid_index.eq.2) hist1_eta_2(index) =
4426 1 hist1_eta_2(index) + 1
4432 if(pid_index .eq. 1) n_part_used_1_trk2 = trk2_counter
4433 if(pid_index .eq. 2) n_part_used_2_trk2 = trk2_counter
4436 End If ! End Mode = 1,2,3,4 Selection Options
4442 C-----------------------------------------------------------------------
4444 subroutine histog2(mode,itrack,px_sec_save,py_sec_save,
4445 1 pz_sec_save,px_save,py_save,pz_save,E_save)
4448 CCC This subroutine computes and/or updates the two-body histograms
4449 C according to the selected 'mode' and for the necessary particle
4452 C Mode = 1, Fill the two-body histograms (hist*) for all particles
4453 C in table 'trk' for like and unlike pairs, for 1D and/or
4454 C 3D fine and 3D coarse mesh bins.
4456 C Mode = 2, Remove all old track pairs for 'itrack' particle from
4457 C all htmp* histograms, for particles in table 'trk', for
4458 C like and unlike pairs, for 1D and/or 3D fine and 3D coarse
4459 C mesh bins; using the saved momentum and track values.
4461 C Mode = 3, Add all new track pairs for 'itrack' particle to
4462 C all htmp* histograms, for particles in table 'trk', for
4463 C like and unlike pairs, for 1D and/or 3D fine and 3D coarse
4464 C mesh bins; using the values in table 'trk(itrack).*'
4466 C Mode = 4, Fill and accumulate reference histograms (href*) for all
4467 C particle pairs from tables 'trk' and 'trk2', for like and
4468 C unlike pairs, for 1D and/or 3D fine and 3D coarse
4471 C itrack = single track index in table 'trk' for pairs to be removed
4472 C (mode = 2) or added (mode = 3).
4474 Include 'common_parameters.inc'
4475 Include 'common_mesh.inc'
4476 Include 'common_histograms.inc'
4478 Include 'common_track.inc'
4479 Include 'common_track2.inc'
4480 Include 'common_sec_track.inc'
4481 Include 'common_sec_track2.inc'
4483 CCC Local Variable Type Declarations:
4485 integer*4 mode,itrack,i,j,k,jx,jy,jz
4486 integer*4 jsec,trkj_sector,imin,imax,njtrks
4487 integer*4 index1,index2,index3
4488 integer*4 findex1,findex2,findex3
4489 integer*4 ipxmin,ipymin,ipzmin
4490 integer*4 ipxmax,ipymax,ipzmax
4491 integer*4 trki_pid,trkj_pid
4492 integer*4 px_sec_save,py_sec_save,pz_sec_save
4494 real*4 qinvar2,qtotal2,qvector2
4495 real*4 qside2, qout2, qlong2
4496 real*4 qperp2, qtime2
4497 real*4 qdiff1, qdiff2, qdiff3
4498 real*4 px_save,py_save,pz_save,E_save
4500 If (mode .eq. 1) Then ! Full hist* filling; initialize arrays to zero.
4504 hist_unlike_1d(i) = 0
4510 hist_like_3d_fine(i,j,k) = 0
4511 hist_unlike_3d_fine(i,j,k) = 0
4512 hist_like_3d_coarse(i,j,k) = 0
4513 hist_unlike_3d_coarse(i,j,k) = 0
4520 CCC Select # of particles to loop over for each 'mode':
4522 If (mode .eq. 1) Then
4524 imax = n_part_tot_trk
4525 Else If (mode .eq. 2 .or. mode .eq. 3) Then
4528 Else If (mode .eq. 4) Then
4530 imax = n_part_tot_trk
4533 C------------------------------------------------------
4534 CCC Begin Primary Loop over particles in Table 'trk':
4535 C------------------------------------------------------
4538 if(trk_flag(i) .eq. 0) then
4539 trki_pid = trk_ge_pid(i)
4541 Call sec_index(px_sec_save,n_px_bins,ipxmin,ipxmax)
4542 Call sec_index(py_sec_save,n_py_bins,ipymin,ipymax)
4543 Call sec_index(pz_sec_save,n_pz_bins,ipzmin,ipzmax)
4545 Call sec_index(trk_px_sec(i),n_px_bins,ipxmin,ipxmax)
4546 Call sec_index(trk_py_sec(i),n_py_bins,ipymin,ipymax)
4547 Call sec_index(trk_pz_sec(i),n_pz_bins,ipzmin,ipzmax)
4550 CCC Begin Loop over neighboring sectors for track # 'i':
4552 do jx = ipxmin,ipxmax
4553 do jy = ipymin,ipymax
4554 do jz = ipzmin,ipzmax
4555 trkj_sector = jx + (jy-1)*n_px_bins
4556 1 + (jz-1)*n_px_bins*n_py_bins
4558 if(mode.le.3) njtrks = stm_n_trk_sec(trkj_sector)
4559 if(mode.eq.4) njtrks = stm2_n_trk_sec(trkj_sector)
4560 if(njtrks .gt. 0) then
4562 CCC Begin Secondary Loop over particles in selected sectors in tables
4563 CCC 'trk' or 'trk2':
4565 if(mode.le.3) j = stm_track_id(jsec,trkj_sector)
4566 if(mode.eq.4) j = stm2_track_id(jsec,trkj_sector)
4567 if((mode.eq.1 .and. j.lt.i .and. trk_flag(j).eq.0)
4568 1 .or.(mode.eq.2 .and. j.ne.i .and. trk_flag(j).eq.0)
4569 2 .or.(mode.eq.3 .and. j.ne.i .and. trk_flag(j).eq.0)
4570 3 .or.(mode.eq.4 .and. trk2_flag(j).eq.0)) then
4572 CCC Obtain 1D and 3D relative momenta:
4574 if(mode.eq.1 .or. mode.eq.3) then
4575 trkj_pid = trk_ge_pid(j)
4576 Call qdiff(trk_px(i),trk_py(i),trk_pz(i),trk_E(i),
4577 1 trk_px(j),trk_py(j),trk_pz(j),trk_E(j),
4578 2 qinvar2,qtotal2,qvector2,qside2,qout2,qlong2,
4580 else if(mode.eq.2) then
4581 trkj_pid = trk_ge_pid(j)
4582 Call qdiff(px_save,py_save,pz_save,E_save,
4583 1 trk_px(j),trk_py(j),trk_pz(j),trk_E(j),
4584 2 qinvar2,qtotal2,qvector2,qside2,qout2,qlong2,
4586 else if(mode.eq.4) then
4587 trkj_pid = trk2_ge_pid(j)
4588 Call qdiff(trk_px(i),trk_py(i),trk_pz(i),trk_E(i),
4589 1 trk2_px(j),trk2_py(j),trk2_pz(j),trk2_E(j),
4590 2 qinvar2,qtotal2,qvector2,qside2,qout2,qlong2,
4594 C-----------------------------------------------
4595 CCC Fill and/or Update 1D two-body Histograms:
4596 C-----------------------------------------------
4598 if(switch_1d .gt. 0) then
4600 if(switch_1d .eq. 1) then
4601 qdiff1 = sqrt(qinvar2)
4602 else if(switch_1d .eq. 2) then
4603 qdiff1 = sqrt(qtotal2)
4604 else if(switch_1d .eq. 3) then
4605 qdiff1 = sqrt(qvector2)
4607 qdiff1 = sqrt(qvector2)
4610 if(qdiff1 .le. qmid_1d) then
4611 index1 = int((qdiff1/binsize_1d_fine)+0.99999)
4612 if(index1 .eq. 0) index1 = 1
4613 else if(qdiff1.gt.qmid_1d.and.qdiff1.le.qmax_1d) then
4614 index1 = int(((qdiff1-qmid_1d)/binsize_1d_coarse)
4616 if(index1.eq.0) index1 = 1
4617 index1 = index1 + n_1d_fine
4622 if(index1.ge.1.and.index1.le.n_1d_total) then
4623 if((trki_pid.eq.trkj_pid).and.(switch_type.eq.1
4624 1 .or. switch_type.eq.3)) then
4626 hist_like_1d(index1) = hist_like_1d(index1) + 1
4627 else if(mode.eq.2) then
4628 htmp_like_1d(index1) = htmp_like_1d(index1) - 1
4629 else if(mode.eq.3) then
4630 htmp_like_1d(index1) = htmp_like_1d(index1) + 1
4631 else if(mode.eq.4) then
4632 href_like_1d(index1) = href_like_1d(index1) + 1
4635 else if((trki_pid.ne.trkj_pid).and.(switch_type.eq.2
4636 1 .or. switch_type.eq.3)) then
4638 hist_unlike_1d(index1) = hist_unlike_1d(index1)+1
4639 else if(mode.eq.2) then
4640 htmp_unlike_1d(index1) = htmp_unlike_1d(index1)-1
4641 else if(mode.eq.3) then
4642 htmp_unlike_1d(index1) = htmp_unlike_1d(index1)+1
4643 else if(mode.eq.4) then
4644 href_unlike_1d(index1) = href_unlike_1d(index1)+1
4649 end if ! End 1D Histogram Fill and/or Update.
4651 C-----------------------------------------------
4652 CCC Fill and/or Update 3D Two-Body Histograms:
4653 C-----------------------------------------------
4655 if(switch_3d .gt. 0) then
4656 if(switch_3d .eq. 1) then
4657 qdiff1 = sqrt(qside2)
4658 qdiff2 = sqrt(qout2)
4659 qdiff3 = sqrt(qlong2)
4660 else if(switch_3d .eq. 2) then
4661 qdiff1 = sqrt(qperp2)
4662 qdiff2 = sqrt(qtime2)
4663 qdiff3 = sqrt(qlong2)
4665 qdiff1 = sqrt(qperp2)
4666 qdiff2 = sqrt(qtime2)
4667 qdiff3 = sqrt(qlong2)
4670 if(qdiff1 .le. qmid_3d) then
4671 findex1 = int((qdiff1/binsize_3d_fine)+0.99999)
4672 if(findex1 .eq. 0) findex1 = 1
4674 else if(qdiff1.gt.qmid_3d.and.qdiff1.le.qmax_3d) then
4675 index1 = int((qdiff1/binsize_3d_coarse)+0.99999)
4676 if(index1.eq.1) index1 = 2
4683 if(qdiff2 .le. qmid_3d) then
4684 findex2 = int((qdiff2/binsize_3d_fine)+0.99999)
4685 if(findex2 .eq. 0) findex2 = 1
4687 else if(qdiff2.gt.qmid_3d.and.qdiff2.le.qmax_3d) then
4688 index2 = int((qdiff2/binsize_3d_coarse)+0.99999)
4689 if(index2.eq.1) index2 = 2
4696 if(qdiff3 .le. qmid_3d) then
4697 findex3 = int((qdiff3/binsize_3d_fine)+0.99999)
4698 if(findex3 .eq. 0) findex3 = 1
4700 else if(qdiff3.gt.qmid_3d.and.qdiff3.le.qmax_3d) then
4701 index3 = int((qdiff3/binsize_3d_coarse)+0.99999)
4702 if(index3.eq.1) index3 = 2
4709 if((index1.ge.1.and.index1.le.n_3d_coarse).and.
4710 1 (index2.ge.1.and.index2.le.n_3d_coarse).and.
4711 2 (index3.ge.1.and.index3.le.n_3d_coarse)) then
4713 if((index1+index2+index3).eq.3) then
4715 if(findex1.ge.1.and.findex1.le.n_3d_fine.and.
4716 1 findex2.ge.1.and.findex2.le.n_3d_fine.and.
4717 2 findex3.ge.1.and.findex3.le.n_3d_fine) then
4719 if((trki_pid.eq.trkj_pid).and.(switch_type.eq.1
4720 1 .or. switch_type.eq.3)) then
4723 hist_like_3d_fine(findex1,findex2,findex3) =
4724 1 hist_like_3d_fine(findex1,findex2,findex3) +1
4725 else if(mode.eq.2) then
4726 htmp_like_3d_fine(findex1,findex2,findex3) =
4727 1 htmp_like_3d_fine(findex1,findex2,findex3) -1
4728 else if(mode.eq.3) then
4729 htmp_like_3d_fine(findex1,findex2,findex3) =
4730 1 htmp_like_3d_fine(findex1,findex2,findex3) +1
4731 else if(mode.eq.4) then
4732 href_like_3d_fine(findex1,findex2,findex3) =
4733 1 href_like_3d_fine(findex1,findex2,findex3) +1
4736 else if((trki_pid.ne.trkj_pid).and.(switch_type
4737 1 .eq.2 .or. switch_type.eq.3)) then
4740 hist_unlike_3d_fine(findex1,findex2,findex3) =
4741 1 hist_unlike_3d_fine(findex1,findex2,findex3) +1
4742 else if(mode.eq.2) then
4743 htmp_unlike_3d_fine(findex1,findex2,findex3) =
4744 1 htmp_unlike_3d_fine(findex1,findex2,findex3) -1
4745 else if(mode.eq.3) then
4746 htmp_unlike_3d_fine(findex1,findex2,findex3) =
4747 1 htmp_unlike_3d_fine(findex1,findex2,findex3) +1
4748 else if(mode.eq.4) then
4749 href_unlike_3d_fine(findex1,findex2,findex3) =
4750 1 href_unlike_3d_fine(findex1,findex2,findex3) +1
4756 else if((index1+index2+index3).gt.3) then
4758 if((trki_pid.eq.trkj_pid).and.(switch_type.eq.1
4759 1 .or. switch_type.eq.3)) then
4762 hist_like_3d_coarse(index1,index2,index3) =
4763 1 hist_like_3d_coarse(index1,index2,index3) +1
4764 else if(mode.eq.2) then
4765 htmp_like_3d_coarse(index1,index2,index3) =
4766 1 htmp_like_3d_coarse(index1,index2,index3) -1
4767 else if(mode.eq.3) then
4768 htmp_like_3d_coarse(index1,index2,index3) =
4769 1 htmp_like_3d_coarse(index1,index2,index3) +1
4770 else if(mode.eq.4) then
4771 href_like_3d_coarse(index1,index2,index3) =
4772 1 href_like_3d_coarse(index1,index2,index3) +1
4775 else if((trki_pid.ne.trkj_pid).and.(switch_type
4776 1 .eq.2 .or. switch_type.eq.3)) then
4779 hist_unlike_3d_coarse(index1,index2,index3) =
4780 1 hist_unlike_3d_coarse(index1,index2,index3) +1
4781 else if(mode.eq.2) then
4782 htmp_unlike_3d_coarse(index1,index2,index3) =
4783 1 htmp_unlike_3d_coarse(index1,index2,index3) -1
4784 else if(mode.eq.3) then
4785 htmp_unlike_3d_coarse(index1,index2,index3) =
4786 1 htmp_unlike_3d_coarse(index1,index2,index3) +1
4787 else if(mode.eq.4) then
4788 href_unlike_3d_coarse(index1,index2,index3) =
4789 1 href_unlike_3d_coarse(index1,index2,index3) +1
4793 end if ! End 3D Fine/Coarse Grid
4795 end if ! End 3D Histogram Filling and/or Update
4798 end do ! End Secondary Track Loop
4803 end do ! End Neighboring Sector Loop
4806 end do ! End Primary Track Loop
4811 C-----------------------------------------------------------------------
4814 subroutine correlation_fit
4817 CCC This subroutine carries out the track momentum adjustment
4818 CCC procedure in order to fit the requested model correlation
4819 CCC function and the input one-body {pt,phi,eta} distributions.
4821 CCC The method used is similar to the Metropolis method. Briefly:
4823 C 1. The accepted tracks for each event in the 'event_text.in'
4824 C input file are loaded into the 'trk' data structure table.
4825 C 2. The sector information, histograms, C2 and initial chi-square
4827 C 3. Each track momentum is randomly shifted within a specified
4828 C range, one track at a time, the histograms are updated, and
4829 C a new C2 and chi-square are computed.
4830 C 4. If the new track momentum is acceptable and if it results in a
4831 C smaller value of chi-square, then this shifted momentum is
4832 C retained, if not then the original momentum value is restored.
4833 C 5. This is done for all particles in the track table for the event.
4834 C 6. The entire process is repeated either until the maximum # of
4835 C iterations is reached, or until the chi-square improvement with
4836 C each iteration diminishes sufficiently.
4837 C 7. After completing the event loop, summary information is gathered
4838 C and inclusive correlation functions and one-body distributions
4841 Include 'common_parameters.inc'
4842 Include 'common_mesh.inc'
4843 Include 'common_histograms.inc'
4844 Include 'common_correlations.inc'
4845 Include 'common_event_summary.inc'
4847 Include 'common_track.inc'
4848 Include 'common_track2.inc'
4849 Include 'common_sec_track.inc'
4850 Include 'common_sec_track2.inc'
4851 Include 'common_particle.inc'
4853 CCC Local Variable Type Declarations:
4855 integer*4 i,j,k,ievent,niter,ntracks_out,nev
4856 integer*4 ntracks_flagged,track_status,pid_index
4857 integer*4 px_sec_save,py_sec_save,pz_sec_save
4859 real*4 px_save,py_save,pz_save,E_save
4860 real*4 pt_save,phi_save,eta_save,mass
4861 real*4 chisq_like_1d,chisq_unlike_1d
4862 real*4 chisq_like_3d_fine,chisq_unlike_3d_fine
4863 real*4 chisq_like_3d_coarse,chisq_unlike_3d_coarse
4864 real*4 chisq_hist1_1,chisq_hist1_2
4865 real*4 chisq_total,chisq_total_oldvalue,chisq_total_newvalue
4868 CCC Initialize counters:
4870 n_part_used_1_inc = 0
4871 n_part_used_2_inc = 0
4872 num_pairs_like_inc = 0
4873 num_pairs_unlike_inc = 0
4874 event_line_counter = 0
4875 file10_line_counter = 0
4877 CCC Open event input, track selection flags and event output files:
4879 If(ALICE .eq. 0) Then
4880 open(unit=2,type='old',access='sequential',
4881 1 name='event_text.in')
4882 open(unit=4,type='old',access='sequential',
4883 1 name='event_tracks.select')
4884 open(unit=10,status='unknown',access='sequential',
4885 1 name='event_hbt_text.out')
4886 CCC Read/Write event header from/to I/O event text files
4891 C-------------------------------------
4894 do ievent = 2, n_events + 1
4895 C-------------------------------------
4897 If(ALICE .eq. 1) Then
4898 Call AliHbtp_SetActiveEventNumber(ievent-1)
4899 C write(*,*) 'NEXT EVENT:', ievent
4902 if(n_part_tot_trk .gt. 0) then
4905 Call tindex(1,0) ! Fill initial track-sector info.
4906 Call stm_build(1,0,0) ! Fill initial sector info.
4907 Call dist_range(1,ntracks_out,ntracks_flagged)
4908 num_pairs_like = (n_part_used_1_trk*(n_part_used_1_trk-1))/2
4909 1 + (n_part_used_2_trk*(n_part_used_2_trk-1))/2
4910 num_pairs_unlike = n_part_used_1_trk*n_part_used_2_trk
4911 num_pairs_like_inc = num_pairs_like_inc + num_pairs_like
4912 num_pairs_unlike_inc = num_pairs_unlike_inc + num_pairs_unlike
4913 n_part_used_1_inc = n_part_used_1_inc + n_part_used_1_trk
4914 n_part_used_2_inc = n_part_used_2_inc + n_part_used_2_trk
4916 C write (*,*) 'num_pairs_like = ',num_pairs_like
4917 C write (*,*) 'num_pairs_unlike = ',num_pairs_unlike
4918 C write (*,*) 'num_pairs_like_inc = ',num_pairs_like_inc
4919 c write (*,*) 'num_pairs_unlike_inc = ',num_pairs_unlike_inc
4920 c write (*,*) 'n_part_used_1_inc = ',n_part_used_1_inc
4921 C write (*,*) 'n_part_used_2_inc = ',n_part_used_2_inc
4923 if(pid(1).gt.0) Call histog1(1,0,1,pid(1),0.,0.,0.)
4924 if(pid(2).gt.0) Call histog1(1,0,2,pid(2),0.,0.,0.)
4925 Call histog2(1,0,0,0,0,0.0,0.0,0.0,0.0)
4927 Call chisquare(1,chisq_like_1d,chisq_unlike_1d,
4928 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
4929 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
4930 3 chisq_hist1_1,chisq_hist1_2)
4931 chisq_total = chisq_wt_like_1d *chisq_like_1d
4932 1 + chisq_wt_unlike_1d *chisq_unlike_1d
4933 2 + chisq_wt_like_3d_fine *chisq_like_3d_fine
4934 3 + chisq_wt_unlike_3d_fine *chisq_unlike_3d_fine
4935 4 + chisq_wt_like_3d_coarse *chisq_like_3d_coarse
4936 5 + chisq_wt_unlike_3d_coarse *chisq_unlike_3d_coarse
4937 6 + chisq_wt_hist1_1 *chisq_hist1_1
4938 7 + chisq_wt_hist1_2 *chisq_hist1_2
4939 chisq_total_oldvalue = chisq_total
4944 1000 Continue ! Starting Point for Track Shift Iteration Loop:
4952 98 Format(5x,'** START NEXT EVENT **')
4953 99 Format(5x,'************************')
4954 write(8,100) ievent,niter,chisq_total
4955 100 Format(10x,'Event#+1,Iteration# and Chi-Sq = ',2I5,E11.4)
4957 IF(maxit .eq. 0) GO TO 1001 ! Option to compute correlations
4958 C ! for input events.
4960 C-------------------------------------
4961 C Begin Track Adjustment Loop:
4963 do i = 1,n_part_tot_trk
4964 C-------------------------------------
4966 if(trk_flag(i) .eq. 0) then
4968 CCC Save initial track parameters (those that might change):
4975 phi_save = trk_phi(i)
4976 eta_save = trk_eta(i)
4977 px_sec_save = trk_px_sec(i)
4978 py_sec_save = trk_py_sec(i)
4979 pz_sec_save = trk_pz_sec(i)
4980 old_sec_save = trk_sector(i)
4982 CCC Save the sector values that the track is initially located in:
4984 old_sec_ntrk = stm_n_trk_sec(trk_sector(i))
4985 old_sec_flag = stm_flag(trk_sector(i))
4986 do k = 1,stm_n_trk_sec(trk_sector(i))
4987 old_sec_trkid(k) = stm_track_id(k,trk_sector(i))
4990 CCC Determine the particle ID type:
4992 if(trk_ge_pid(i).eq.pid(1) .and. pid(1).gt.0) then
4994 else if(trk_ge_pid(i).eq.pid(2).and.pid(2).gt.0) then
5000 CCC Randomly shift track momentum vector and compute new kinematics:
5002 trk_px(i) = trk_px(i) + deltap*(2.0*hbtpran(irand) - 1.0)
5003 trk_py(i) = trk_py(i) + deltap*(2.0*hbtpran(irand) - 1.0)
5004 trk_pz(i) = trk_pz(i) + deltap*(2.0*hbtpran(irand) - 1.0)
5005 mass = part_mass(trk_ge_pid(i))
5006 Call Hbtp_kin(trk_px(i),trk_py(i),trk_pz(i),trk_E(i),
5007 1 trk_pt(i),trk_phi(i),trk_eta(i),mass,1)
5009 new_sec_save = trk_sector(i)
5011 CCC Determine if track has been shifted to a new sector, and if so,
5012 CCC whether this overfills this new sector. If all is well, then
5013 CCC update histograms. If not, then restore track parameters to their
5014 CCC initial values prior to shifting. Keep the status of track(i) in
5015 CCC 'track_status', where a value of 0 means the track is OK to use.
5017 CCC The Logical steps are the following:
5019 C IF(new track position is in same sector) THEN
5020 C o Remove old track position from htmp1*, htmp* using old saved values.
5021 C o Add new track position to htmp1*, htmp* using values in 'trk'
5022 C (Sector information is unchanged)
5023 C ELSE IF(new track position is in a different sector) THEN
5024 C IF(# tracks in new sector is still OK, with the new track) THEN
5025 C o Save values of new sector before trk was shifted into it.
5026 C o Remove old trk position from htmp1*, htmp* using old saved values
5027 C o Add new trk position to htmp1*, htmp* using values in trk
5028 C o Update sector information in stm
5029 C ELSE IF(# tracks in new sector becomes too many with new trk) THEN
5030 C o Restore track parameters to pre-shifted values
5031 C o Set track_status = 1, indicating the track could not be moved
5036 if(old_sec_save .eq. new_sec_save) then
5037 Call histog1(2,i,pid_index,pid(pid_index),pt_save,
5038 1 phi_save,eta_save)
5039 Call histog2(2,i,px_sec_save,py_sec_save,pz_sec_save,
5040 1 px_save,py_save,pz_save,E_save)
5042 Call histog1(3,i,pid_index,pid(pid_index),0.,0.,0.)
5043 Call histog2(3,i,0,0,0,0.0,0.0,0.0,0.0)
5045 else if(old_sec_save .ne. new_sec_save) then
5047 if(stm_n_trk_sec(new_sec_save) .lt. max_trk_sec) then
5048 new_sec_ntrk = stm_n_trk_sec(new_sec_save)
5049 new_sec_flag = stm_flag(new_sec_save)
5050 if(new_sec_ntrk .gt. 0) then
5051 do k = 1,new_sec_ntrk
5052 new_sec_trkid(k) = stm_track_id(k,new_sec_save)
5056 Call histog1(2,i,pid_index,pid(pid_index),pt_save,
5057 1 phi_save,eta_save)
5058 Call histog2(2,i,px_sec_save,py_sec_save,pz_sec_save,
5059 1 px_save,py_save,pz_save,E_save)
5061 Call histog1(3,i,pid_index,pid(pid_index),
5063 Call histog2(3,i,0,0,0,0.0,0.0,0.0,0.0)
5065 Call stm_build(1,i,old_sec_save)
5067 else if(stm_n_trk_sec(new_sec_save) .ge. max_trk_sec) then
5075 trk_phi(i) = phi_save
5076 trk_eta(i) = eta_save
5077 trk_px_sec(i) = px_sec_save
5078 trk_py_sec(i) = py_sec_save
5079 trk_pz_sec(i) = pz_sec_save
5080 trk_sector(i) = old_sec_save
5083 end if ! End Histogram and Sector Update
5085 CCC If the track was succesfully shifted then compute C2 and determine
5086 C if the chi-square decreases (improves) or increases. If it improves,
5087 C then store the new chi-square value and keep the 1- and 2-body
5088 C histograms in hist1* and hist*, repsectively. If chi-square
5089 C increases (worsens), then restore the track parameters to the
5090 C pre-shifted values, restore the histograms and if a new sector was
5091 C involved, then restore both the old and new sector values.
5093 C The Logical steps are the following:
5095 C IF(new track position is OK, (i.e. track_status = 0)) Then
5096 C o Compute C2 using htmp*
5097 C o Compute chi-square and save
5098 C IF(chi-square improves) Then
5099 C o Replace previous (best) chi-square with new value
5100 C o Update histograms, i.e. copy htmp1* -> hist1* and
5101 C copy htmp* -> hist*
5102 C ELSE IF(chi-square increases) Then
5103 C o Restore track parameters
5104 C o Restore histograms, i.e. copy hist1* -> htmp1* and
5105 C copy hist* -> htmp*
5106 C IF(new sector was used) Then
5107 C o Restore old sector values to pre-shifted values
5108 C o Restore new sector values to pre-shifted values
5113 If(track_status .eq.0) Then
5115 Call chisquare(2,chisq_like_1d,chisq_unlike_1d,
5116 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
5117 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
5118 3 chisq_hist1_1,chisq_hist1_2)
5119 chisq_total_newvalue =
5120 1 chisq_wt_like_1d *chisq_like_1d
5121 2 + chisq_wt_unlike_1d *chisq_unlike_1d
5122 3 + chisq_wt_like_3d_fine *chisq_like_3d_fine
5123 4 + chisq_wt_unlike_3d_fine *chisq_unlike_3d_fine
5124 5 + chisq_wt_like_3d_coarse *chisq_like_3d_coarse
5125 6 + chisq_wt_unlike_3d_coarse *chisq_unlike_3d_coarse
5126 7 + chisq_wt_hist1_1 *chisq_hist1_1
5127 8 + chisq_wt_hist1_2 *chisq_hist1_2
5129 if(chisq_total_newvalue .lt. chisq_total_oldvalue) then
5130 chisq_total_oldvalue = chisq_total_newvalue
5133 else if(chisq_total_newvalue.ge.chisq_total_oldvalue) then
5139 trk_phi(i) = phi_save
5140 trk_eta(i) = eta_save
5141 trk_px_sec(i) = px_sec_save
5142 trk_py_sec(i) = py_sec_save
5143 trk_pz_sec(i) = pz_sec_save
5144 trk_sector(i) = old_sec_save
5148 If(old_sec_save .ne. new_sec_save) then
5150 stm_n_trk_sec(old_sec_save) = old_sec_ntrk
5151 stm_flag(old_sec_save) = old_sec_flag
5152 do k = 1,max_trk_sec
5153 stm_track_id(k,old_sec_save) = 0
5155 do k = 1,old_sec_ntrk
5156 stm_track_id(k,old_sec_save) = old_sec_trkid(k)
5159 stm_n_trk_sec(new_sec_save) = new_sec_ntrk
5160 stm_flag(new_sec_save) = new_sec_flag
5161 do k = 1,max_trk_sec
5162 stm_track_id(k,new_sec_save) = 0
5164 do k = 1,new_sec_ntrk
5165 stm_track_id(k,new_sec_save) = new_sec_trkid(k)
5170 end if ! End Chi-Square Determination
5171 end if ! End valid tracks option
5172 end do ! End Track Shift Loop
5174 CCC Check chi-square for this iteration --
5175 C Best, current chi-square value is in 'chisq_total_oldvalue'
5176 C Chi-square value at the beginning of the iteration loop is in
5179 If(abs(200.0*(chisq_total_oldvalue - chisq_total)/
5180 1 (chisq_total_oldvalue + chisq_total)) .lt. delchi) then
5182 101 Format(/5x,'Chi-Sq reduced .lt. delchi % on last iteration',
5186 If (niter .gt. maxit) Then
5188 102 Format(/5x,'Max # Search Iterations Reached - Abort track ',
5192 chisq_total = chisq_total_oldvalue
5197 CCC Finished Track Adjustment Iteration Loop for event # 'ievent'
5199 if((ievent - 1) .le. max_events) then
5200 Call dist_range(1,ntracks_out,ntracks_flagged)
5201 num_iter(ievent-1) = float(niter)
5202 n_part_used_1_store(ievent-1) = float(n_part_used_1_trk)
5203 n_part_used_2_store(ievent-1) = float(n_part_used_2_trk)
5204 n_part_tot_store(ievent-1) = float(n_part_tot_trk)
5205 frac_trks_out(ievent-1)=float(ntracks_out)/
5206 1 float(n_part_tot_trk)
5207 frac_trks_flag(ievent-1) =
5208 1 float(ntracks_flagged)/float(n_part_tot_trk)
5212 Call chisquare(1,chisq_like_1d,chisq_unlike_1d,
5213 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
5214 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
5215 3 chisq_hist1_1,chisq_hist1_2)
5216 chisq_total = chisq_wt_like_1d *chisq_like_1d
5217 1 + chisq_wt_unlike_1d *chisq_unlike_1d
5218 2 + chisq_wt_like_3d_fine *chisq_like_3d_fine
5219 3 + chisq_wt_unlike_3d_fine *chisq_unlike_3d_fine
5220 4 + chisq_wt_like_3d_coarse *chisq_like_3d_coarse
5221 5 + chisq_wt_unlike_3d_coarse *chisq_unlike_3d_coarse
5222 6 + chisq_wt_hist1_1 *chisq_hist1_1
5223 7 + chisq_wt_hist1_2 *chisq_hist1_2
5225 if((ievent - 1) .le. max_events) then
5226 chisq_like_1d_store(ievent-1) = chisq_like_1d
5227 chisq_unlike_1d_store(ievent-1) = chisq_unlike_1d
5228 chisq_like_3d_fine_store(ievent-1) = chisq_like_3d_fine
5229 chisq_unlike_3d_fine_store(ievent-1) = chisq_unlike_3d_fine
5230 chisq_like_3d_coarse_store(ievent-1) = chisq_like_3d_coarse
5231 chisq_unlike_3d_coarse_store(ievent-1) =chisq_unlike_3d_coarse
5232 chisq_hist1_1_store(ievent-1) = chisq_hist1_1
5233 chisq_hist1_2_store(ievent-1) = chisq_hist1_2
5234 chisq_total_store(ievent-1) = chisq_total
5236 CCC Count # sectors with stm().flag = 1, indicating that too many
5237 C tracks were attempted to be loaded into that sector.
5239 num_sec_flagged_store(ievent-1) = 0.0
5241 if(stm_flag(k) .eq. 1) then
5242 num_sec_flagged_store(ievent-1) =
5243 1 num_sec_flagged_store(ievent-1) + 1.0
5250 if(print_full .eq. 1) Call write_data(5,ievent-1)
5252 end if ! End event-with-tracks processing.
5255 C-------------------------------
5256 end do ! End Event Loop
5257 C-------------------------------
5259 CCC Compute Correlation Functions for the Inclusive Histograms
5263 CCC Compute Mean and Std. dev of event monitor and summary quantities:
5265 if(n_events .le. max_events) then
5271 Call mean_rms(num_iter,nev,nev,niter_mean,niter_rms)
5272 Call mean_rms(n_part_used_1_store,nev,nev,npart1_mean,npart1_rms)
5273 Call mean_rms(n_part_used_2_store,nev,nev,npart2_mean,npart2_rms)
5274 Call mean_rms(n_part_tot_store,nev,nev,npart_tot_mean,
5276 Call mean_rms(num_sec_flagged_store,nev,nev,
5277 1 nsec_flag_mean,nsec_flag_rms)
5278 Call mean_rms(frac_trks_out,nev,nev,
5279 1 frac_trks_out_mean,frac_trks_out_rms)
5280 Call mean_rms(frac_trks_flag,nev,nev,
5281 1 frac_trks_flag_mean,frac_trks_flag_rms)
5282 Call mean_rms(chisq_like_1d_store,nev,nev,
5283 1 chi_l1d_mean,chi_l1d_rms)
5284 Call mean_rms(chisq_unlike_1d_store,nev,nev,
5285 1 chi_u1d_mean,chi_u1d_rms)
5286 Call mean_rms(chisq_like_3d_fine_store,nev,nev,
5287 1 chi_l3f_mean,chi_l3f_rms)
5288 Call mean_rms(chisq_unlike_3d_fine_store,nev,nev,
5289 1 chi_u3f_mean,chi_u3f_rms)
5290 Call mean_rms(chisq_like_3d_coarse_store,nev,nev,
5291 1 chi_l3c_mean,chi_l3c_rms)
5292 Call mean_rms(chisq_unlike_3d_coarse_store,nev,nev,
5293 1 chi_u3c_mean,chi_u3c_rms)
5294 Call mean_rms(chisq_hist1_1_store,nev,nev,
5295 1 chi_1_1_mean, chi_1_1_rms)
5296 Call mean_rms(chisq_hist1_2_store,nev,nev,
5297 1 chi_1_2_mean, chi_1_2_rms)
5298 Call mean_rms(chisq_total_store,nev,nev,
5299 1 chi_tot_mean, chi_tot_rms)
5301 If(ALICE .eq. 0) Then
5310 C------------------------------------------------------------------------
5313 subroutine hist1_copy(mode)
5316 CCC Copy 1-body histograms if:
5318 CCC mode = 1, then copy hist1* -> htmp1*
5319 CCC mode = 2, then copy htmp1* -> hist1*
5321 Include 'common_parameters.inc'
5322 Include 'common_mesh.inc'
5323 Include 'common_histograms.inc'
5325 CCC Local Variable Type Declarations:
5329 C---------------------------
5330 If(mode .eq. 1) Then ! Copy hist1* -> htmp1*
5331 C---------------------------
5333 if(pid(1) .gt. 0) then
5335 htmp1_pt_1(i) = hist1_pt_1(i)
5338 htmp1_phi_1(i) = hist1_phi_1(i)
5341 htmp1_eta_1(i) = hist1_eta_1(i)
5345 if(pid(2) .gt. 0) then
5347 htmp1_pt_2(i) = hist1_pt_2(i)
5350 htmp1_phi_2(i) = hist1_phi_2(i)
5353 htmp1_eta_2(i) = hist1_eta_2(i)
5357 C--------------------------------
5358 Else If (mode .eq. 2) Then ! Copy htmp1* -> hist1*
5359 C--------------------------------
5361 if(pid(1) .gt. 0) then
5363 hist1_pt_1(i) = htmp1_pt_1(i)
5366 hist1_phi_1(i) = htmp1_phi_1(i)
5369 hist1_eta_1(i) = htmp1_eta_1(i)
5373 if(pid(2) .gt. 0) then
5375 hist1_pt_2(i) = htmp1_pt_2(i)
5378 hist1_phi_2(i) = htmp1_phi_2(i)
5381 hist1_eta_2(i) = htmp1_eta_2(i)
5392 C----------------------------------------------------------------------
5394 subroutine hist2_copy(mode)
5397 CCC Copy 2-body histograms if:
5399 CCC mode = 1, then copy hist* -> htmp*
5400 CCC mode = 2, then copy htmp* -> hist*
5402 Include 'common_parameters.inc'
5403 Include 'common_mesh.inc'
5404 Include 'common_histograms.inc'
5406 CCC Local Variable Type Declarations:
5408 integer*4 mode, i,j,k
5410 C---------------------------
5411 If (mode .eq. 1) Then ! Copy hist* -> htmp*
5412 C---------------------------
5414 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
5415 if(switch_type.eq.1 .or. switch_type.eq.3) then
5417 htmp_like_1d(i) = hist_like_1d(i)
5420 if(switch_type.eq.2 .or. switch_type.eq.3) then
5422 htmp_unlike_1d(i) = hist_unlike_1d(i)
5425 end if ! End 1D histogram copy
5427 if(switch_3d.gt.0) then
5428 if(switch_type.eq.1 .or. switch_type.eq.3) then
5430 if(n_3d_fine .gt. 0) then
5434 htmp_like_3d_fine(i,j,k) = hist_like_3d_fine(i,j,k)
5440 if(n_3d_coarse .gt. 0) then
5441 do i = 1,n_3d_coarse
5442 do j = 1,n_3d_coarse
5443 do k = 1,n_3d_coarse
5444 htmp_like_3d_coarse(i,j,k) = hist_like_3d_coarse(i,j,k)
5452 if(switch_type.eq.2 .or. switch_type.eq.3) then
5454 if(n_3d_fine .gt. 0) then
5458 htmp_unlike_3d_fine(i,j,k) = hist_unlike_3d_fine(i,j,k)
5464 if(n_3d_coarse .gt. 0) then
5465 do i = 1,n_3d_coarse
5466 do j = 1,n_3d_coarse
5467 do k = 1,n_3d_coarse
5468 htmp_unlike_3d_coarse(i,j,k)=hist_unlike_3d_coarse(i,j,k)
5475 end if ! End 3D histogram copy
5477 C--------------------------------
5478 Else If (mode .eq. 2) Then ! Copy htmp* -> hist*
5479 C--------------------------------
5481 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
5482 if(switch_type.eq.1 .or. switch_type.eq.3) then
5484 hist_like_1d(i) = htmp_like_1d(i)
5487 if(switch_type.eq.2 .or. switch_type.eq.3) then
5489 hist_unlike_1d(i) = htmp_unlike_1d(i)
5492 end if ! End 1D histogram copy
5494 if(switch_3d.gt.0) then
5495 if(switch_type.eq.1 .or. switch_type.eq.3) then
5497 if(n_3d_fine .gt. 0) then
5501 hist_like_3d_fine(i,j,k) = htmp_like_3d_fine(i,j,k)
5507 if(n_3d_coarse .gt. 0) then
5508 do i = 1,n_3d_coarse
5509 do j = 1,n_3d_coarse
5510 do k = 1,n_3d_coarse
5511 hist_like_3d_coarse(i,j,k) = htmp_like_3d_coarse(i,j,k)
5519 if(switch_type.eq.2 .or. switch_type.eq.3) then
5521 if(n_3d_fine .gt. 0) then
5525 hist_unlike_3d_fine(i,j,k) = htmp_unlike_3d_fine(i,j,k)
5531 if(n_3d_coarse .gt. 0) then
5532 do i = 1,n_3d_coarse
5533 do j = 1,n_3d_coarse
5534 do k = 1,n_3d_coarse
5535 hist_unlike_3d_coarse(i,j,k)=htmp_unlike_3d_coarse(i,j,k)
5542 end if ! End 3D histogram copy
5545 End If ! End mode selection options
5551 C-----------------------------------------------------------------------
5554 subroutine hist1_incl_sum
5557 CCC Sum 1-body histograms for each event into inclusive totals, where
5558 CCC hinc1* = SUM[hist1*]
5560 Include 'common_parameters.inc'
5561 Include 'common_mesh.inc'
5562 Include 'common_histograms.inc'
5564 CCC Local Variable Type Declarations:
5568 if(pid(1) .gt. 0) then
5570 hinc1_pt_1(i) = hinc1_pt_1(i) + hist1_pt_1(i)
5573 hinc1_phi_1(i) = hinc1_phi_1(i) + hist1_phi_1(i)
5576 hinc1_eta_1(i) = hinc1_eta_1(i) + hist1_eta_1(i)
5580 if(pid(2) .gt. 0) then
5582 hinc1_pt_2(i) = hinc1_pt_2(i) + hist1_pt_2(i)
5585 hinc1_phi_2(i) = hinc1_phi_2(i) + hist1_phi_2(i)
5588 hinc1_eta_2(i) = hinc1_eta_2(i) + hist1_eta_2(i)
5596 C------------------------------------------------------------------------
5599 subroutine hist2_incl_sum
5602 CCC Sum 2-body histograms for each event into inclusive totals, where
5603 CCC hinc* = SUM[hist*]
5605 Include 'common_parameters.inc'
5606 Include 'common_mesh.inc'
5607 Include 'common_histograms.inc'
5609 CCC Local Variable Type Declarations:
5613 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
5614 if(switch_type.eq.1 .or. switch_type.eq.3) then
5616 hinc_like_1d(i) = hinc_like_1d(i) + hist_like_1d(i)
5619 if(switch_type.eq.2 .or. switch_type.eq.3) then
5621 hinc_unlike_1d(i) = hinc_unlike_1d(i) + hist_unlike_1d(i)
5624 end if ! End 1D Inclusive Histogram Sum
5626 if(switch_3d.gt.0) then
5627 if(switch_type.eq.1 .or. switch_type.eq.3) then
5629 if(n_3d_fine .gt. 0) then
5633 hinc_like_3d_fine(i,j,k) = hinc_like_3d_fine(i,j,k)
5634 1 + hist_like_3d_fine(i,j,k)
5640 if(n_3d_coarse .gt. 0) then
5641 do i = 1,n_3d_coarse
5642 do j = 1,n_3d_coarse
5643 do k = 1,n_3d_coarse
5644 hinc_like_3d_coarse(i,j,k) = hinc_like_3d_coarse(i,j,k)
5645 1 + hist_like_3d_coarse(i,j,k)
5653 if(switch_type.eq.2 .or. switch_type.eq.3) then
5655 if(n_3d_fine .gt. 0) then
5659 hinc_unlike_3d_fine(i,j,k) = hinc_unlike_3d_fine(i,j,k)
5660 1 + hist_unlike_3d_fine(i,j,k)
5666 if(n_3d_coarse .gt. 0) then
5667 do i = 1,n_3d_coarse
5668 do j = 1,n_3d_coarse
5669 do k = 1,n_3d_coarse
5670 hinc_unlike_3d_coarse(i,j,k) = hinc_unlike_3d_coarse(i,j,k)
5671 1 + hist_unlike_3d_coarse(i,j,k)
5678 end if ! End 3D Inclusive Histogram Sum
5683 C--------------------------------------------------------------------------
5686 subroutine correl_fit(mode)
5689 CCC This subroutine calculates the 2-body correlation function with
5690 CCC errors for the cases:
5692 CCC (1) 1D and/or 3D fine and coarse grid distributions
5693 CCC (2) like pairs and/or unlike pairs
5695 CCC Uses the signal and reference histograms. The input parameter
5696 CCC 'mode' selects which histograms to use.
5698 CCC Mode = 1, use hist*
5699 CCC Mode = 2, use htmp*
5700 CCC Mode = 3, use hinc*
5702 Include 'common_parameters.inc'
5703 Include 'common_mesh.inc'
5704 Include 'common_histograms.inc'
5705 Include 'common_correlations.inc'
5707 CCC Local Variable Type Declarations:
5709 integer*4 mode,i,j,k
5711 CCC Initialize correlation functions and error arrays to zero:
5714 c2fit_like_1d(i) = 0.0
5715 c2fit_unlike_1d(i) = 0.0
5716 c2err_like_1d(i) = 0.0
5717 c2err_unlike_1d(i) = 0.0
5723 c2fit_like_3d_fine(i,j,k) = 0.0
5724 c2fit_unlike_3d_fine(i,j,k) = 0.0
5725 c2fit_like_3d_coarse(i,j,k) = 0.0
5726 c2fit_unlike_3d_coarse(i,j,k) = 0.0
5727 c2err_like_3d_fine(i,j,k) = 0.0
5728 c2err_unlike_3d_fine(i,j,k) = 0.0
5729 c2err_like_3d_coarse(i,j,k) = 0.0
5730 c2err_unlike_3d_coarse(i,j,k) = 0.0
5735 CCC Compute 1D Correlation Functions and Errors:
5737 if(switch_1d .gt. 0) then
5738 if(switch_type.eq.1 .or. switch_type.eq.3) then
5740 if(mode .eq. 1) then
5741 Call c2_1d(hist_like_1d,href_like_1d,c2fit_like_1d,
5742 1 c2err_like_1d,max_h_1d,max_c2_1d,n_1d_total,
5743 2 num_pairs_like,num_pairs_like_ref)
5744 else if (mode .eq. 2) then
5745 Call c2_1d(htmp_like_1d,href_like_1d,c2fit_like_1d,
5746 1 c2err_like_1d,max_h_1d,max_c2_1d,n_1d_total,
5747 2 num_pairs_like,num_pairs_like_ref)
5748 else if (mode .eq. 3) then
5749 Call c2_1d(hinc_like_1d,href_like_1d,c2fit_like_1d,
5750 1 c2err_like_1d,max_h_1d,max_c2_1d,n_1d_total,
5751 2 num_pairs_like_inc,num_pairs_like_ref)
5756 if(switch_type.eq.2 .or. switch_type.eq.3) then
5758 if(mode .eq. 1) then
5759 Call c2_1d(hist_unlike_1d,href_unlike_1d,c2fit_unlike_1d,
5760 1 c2err_unlike_1d,max_h_1d,max_c2_1d,n_1d_total,
5761 2 num_pairs_unlike,num_pairs_unlike_ref)
5762 else if (mode .eq. 2) then
5763 Call c2_1d(htmp_unlike_1d,href_unlike_1d,c2fit_unlike_1d,
5764 1 c2err_unlike_1d,max_h_1d,max_c2_1d,n_1d_total,
5765 2 num_pairs_unlike,num_pairs_unlike_ref)
5766 else if (mode .eq. 3) then
5767 Call c2_1d(hinc_unlike_1d,href_unlike_1d,c2fit_unlike_1d,
5768 1 c2err_unlike_1d,max_h_1d,max_c2_1d,n_1d_total,
5769 2 num_pairs_unlike_inc,num_pairs_unlike_ref)
5772 end if ! End 1D correlations
5774 CCC Compute 3D Correlation Functions and Errors:
5776 if(switch_3d .gt. 0) then
5777 if(switch_type.eq.1 .or. switch_type.eq.3) then
5779 if(mode .eq. 1) then
5780 Call c2_3d(hist_like_3d_fine,href_like_3d_fine,
5781 1 c2fit_like_3d_fine,c2err_like_3d_fine,
5782 2 max_h_3d,max_c2_3d,n_3d_fine,
5783 3 num_pairs_like,num_pairs_like_ref)
5784 Call c2_3d(hist_like_3d_coarse,href_like_3d_coarse,
5785 1 c2fit_like_3d_coarse,c2err_like_3d_coarse,
5786 2 max_h_3d,max_c2_3d,n_3d_coarse,
5787 3 num_pairs_like,num_pairs_like_ref)
5788 else if(mode .eq. 2) then
5789 Call c2_3d(htmp_like_3d_fine,href_like_3d_fine,
5790 1 c2fit_like_3d_fine,c2err_like_3d_fine,
5791 2 max_h_3d,max_c2_3d,n_3d_fine,
5792 3 num_pairs_like,num_pairs_like_ref)
5793 Call c2_3d(htmp_like_3d_coarse,href_like_3d_coarse,
5794 1 c2fit_like_3d_coarse,c2err_like_3d_coarse,
5795 2 max_h_3d,max_c2_3d,n_3d_coarse,
5796 3 num_pairs_like,num_pairs_like_ref)
5797 else if(mode .eq. 3) then
5798 Call c2_3d(hinc_like_3d_fine,href_like_3d_fine,
5799 1 c2fit_like_3d_fine,c2err_like_3d_fine,
5800 2 max_h_3d,max_c2_3d,n_3d_fine,
5801 3 num_pairs_like_inc,num_pairs_like_ref)
5802 Call c2_3d(hinc_like_3d_coarse,href_like_3d_coarse,
5803 1 c2fit_like_3d_coarse,c2err_like_3d_coarse,
5804 2 max_h_3d,max_c2_3d,n_3d_coarse,
5805 3 num_pairs_like_inc,num_pairs_like_ref)
5810 if(switch_type.eq.2 .or. switch_type.eq.3) then
5812 if(mode .eq. 1) then
5813 Call c2_3d(hist_unlike_3d_fine,href_unlike_3d_fine,
5814 1 c2fit_unlike_3d_fine,c2err_unlike_3d_fine,
5815 2 max_h_3d,max_c2_3d,n_3d_fine,
5816 3 num_pairs_unlike,num_pairs_unlike_ref)
5817 Call c2_3d(hist_unlike_3d_coarse,href_unlike_3d_coarse,
5818 1 c2fit_unlike_3d_coarse,c2err_unlike_3d_coarse,
5819 2 max_h_3d,max_c2_3d,n_3d_coarse,
5820 3 num_pairs_unlike,num_pairs_unlike_ref)
5821 else if(mode .eq. 2) then
5822 Call c2_3d(htmp_unlike_3d_fine,href_unlike_3d_fine,
5823 1 c2fit_unlike_3d_fine,c2err_unlike_3d_fine,
5824 2 max_h_3d,max_c2_3d,n_3d_fine,
5825 3 num_pairs_unlike,num_pairs_unlike_ref)
5826 Call c2_3d(htmp_unlike_3d_coarse,href_unlike_3d_coarse,
5827 1 c2fit_unlike_3d_coarse,c2err_unlike_3d_coarse,
5828 2 max_h_3d,max_c2_3d,n_3d_coarse,
5829 3 num_pairs_unlike,num_pairs_unlike_ref)
5830 else if(mode .eq. 3) then
5831 Call c2_3d(hinc_unlike_3d_fine,href_unlike_3d_fine,
5832 1 c2fit_unlike_3d_fine,c2err_unlike_3d_fine,
5833 2 max_h_3d,max_c2_3d,n_3d_fine,
5834 3 num_pairs_unlike_inc,num_pairs_unlike_ref)
5835 Call c2_3d(hinc_unlike_3d_coarse,href_unlike_3d_coarse,
5836 1 c2fit_unlike_3d_coarse,c2err_unlike_3d_coarse,
5837 2 max_h_3d,max_c2_3d,n_3d_coarse,
5838 3 num_pairs_unlike_inc,num_pairs_unlike_ref)
5841 end if ! End 3D correlations
5847 C-----------------------------------------------------------------------
5850 subroutine c2_1d(h,href,c2,c2err,maxh,maxc2,n,num_pairs_sig,
5854 CCC Computes the two-body correlation function for 1D distributions.
5855 CCC Errors are also computed.
5857 CCC Description of Input Variables in Argument List:
5859 C h(maxh) = signal histogram (numerator)
5860 C href(maxh) = background histogram (denominator)
5861 C c2(maxc2) = correlation function = a/b
5862 C c2err(maxc2) = correlation function error
5863 C maxh = dimension of histogram arrays
5864 C maxc2 = dimension of correlation function array.
5866 C num_pairs_sig = # pairs used in signal histogram
5867 C num_pairs_bkg = # pairs used in background histogram
5870 CCC Local Variable Type Declarations:
5872 integer*4 maxh,maxc2,n,num_pairs_sig,num_pairs_bkg
5873 integer*4 h(maxh), href(maxh)
5876 real*4 c2(maxc2), c2err(maxc2)
5877 real*4 a,a_error,b,b_error
5880 if(href(k).le.0 .or. h(k).le.0) then
5884 a = float(h(k))/float(num_pairs_sig)
5885 a_error = sqrt(float(h(k)))/float(num_pairs_sig)
5886 b = float(href(k))/float(num_pairs_bkg)
5887 b_error = sqrt(float(href(k)))/float(num_pairs_bkg)
5889 c2err(k) = c2(k)*sqrt((a_error/a)**2 + (b_error/b)**2)
5896 C-----------------------------------------------------------------------
5899 subroutine c2_3d(h,href,c2,c2err,maxh,maxc2,n,num_pairs_sig,
5903 CCC Computes the two-body correlation function for 3D distributions.
5904 CCC Errors are also computed.
5906 CCC Description of Input Variables in Argument List:
5908 C h(maxh,maxh,maxh) = 3D signal histogram (numerator)
5909 C href(maxh,maxh,maxh)) = 3D background histogram (denominator)
5910 C c2(maxc2,maxc2,maxc2) = 3D correlation function = a/b
5911 C c2err(maxc2,maxc2,maxc2) = 3D correlation function error
5912 C maxh = dimension of 3D histogram arrays
5913 C maxc2 = dimension of 3D correlation function array.
5915 C num_pairs_sig = # pairs used in signal histogram
5916 C num_pairs_bkg = # pairs used in background histogram
5919 CCC Local Variable Type Declarations:
5921 integer*4 maxh,maxc2,n,num_pairs_sig,num_pairs_bkg
5922 integer*4 h(maxh,maxh,maxh), href(maxh,maxh,maxh)
5925 real*4 c2(maxc2,maxc2,maxc2), c2err(maxc2,maxc2,maxc2)
5926 real*4 a,a_error,b,b_error
5931 if(href(i,j,k).le.0 .or. h(i,j,k).le.0) then
5935 a = float(h(i,j,k))/float(num_pairs_sig)
5936 a_error = sqrt(float(h(i,j,k)))/float(num_pairs_sig)
5937 b = float(href(i,j,k))/float(num_pairs_bkg)
5938 b_error = sqrt(float(href(i,j,k)))/float(num_pairs_bkg)
5940 c2err(i,j,k) = c2(i,j,k)*sqrt((a_error/a)**2 + (b_error/b)**2)
5950 C-------------------------------------------------------------------------
5953 subroutine chisquare(mode,chisq_like_1d,chisq_unlike_1d,
5954 1 chisq_like_3d_fine,chisq_unlike_3d_fine,
5955 2 chisq_like_3d_coarse,chisq_unlike_3d_coarse,
5956 3 chisq_hist1_1,chisq_hist1_2)
5959 CCC This subroutine calculates the chi-squares for the following:
5960 C o Like pair 1D 2-body correlations
5961 C o Unlike pair 1D 2-body correlations
5962 C o Like pair 3D, Fine Mesh 2-body correlations
5963 C o Unlike pair 3D, Fine Mesh 2-body correlations
5964 C o Like pair 3D, Coarse Mesh 2-body correlations
5965 C o Unlike pair 3D, Coarse Mesh 2-body correlations
5966 C o One-body 1D {pt,phi,eta} (summed) distributions for PID#1
5967 C o One-body 1D {pt,phi,eta} (summed) distributions for PID#2
5969 C (where the separate chi-squares for the 1D pt, phi and eta
5970 C one-body distributions are added and only the sum is returned.)
5972 C 'Mode' determines which one-body histogram is compared to the
5973 C reference histogram, where:
5975 C If mode = 1, then hist1* are used
5976 C If mode = 2, then htmp1* are used
5978 C The one-body reference histograms used in the chi-square calculation
5981 Include 'common_parameters.inc'
5982 Include 'common_mesh.inc'
5983 Include 'common_histograms.inc'
5984 Include 'common_correlations.inc'
5986 CCC Local Variable Type Declarations:
5988 integer*4 i,j,k,mode
5990 real*4 chisq_like_1d, chisq_unlike_1d
5991 real*4 chisq_like_3d_fine,chisq_unlike_3d_fine
5992 real*4 chisq_like_3d_coarse,chisq_unlike_3d_coarse
5993 real*4 chisq_hist1_1,chisq_hist1_2
5995 real*4 n1fac,n2fac ! # part 1(2) used/# part 1(2) used in Ref.
5996 real*4 avgerrsq_pt_1, avgerrsq_phi_1, avgerrsq_eta_1
5997 real*4 avgerrsq_pt_2, avgerrsq_phi_2, avgerrsq_eta_2
5998 real*4 avgerrsq_pt_ref_1,avgerrsq_phi_ref_1,avgerrsq_eta_ref_1
5999 real*4 avgerrsq_pt_ref_2,avgerrsq_phi_ref_2,avgerrsq_eta_ref_2
6002 CCC Initialize all chi-square values to zero:
6005 chisq_unlike_1d = 0.0
6006 chisq_like_3d_fine = 0.0
6007 chisq_unlike_3d_fine = 0.0
6008 chisq_like_3d_coarse = 0.0
6009 chisq_unlike_3d_coarse = 0.0
6013 if(switch_1d .gt. 0) then
6014 if(switch_type.eq.1 .or. switch_type.eq.3) then
6016 if(c2fit_like_1d(i) .ne. 0.0) then
6017 chisq_like_1d = chisq_like_1d + ((c2fit_like_1d(i)
6018 1 - c2mod_like_1d(i))/c2err_like_1d(i))**2
6022 if(switch_type.eq.2 .or. switch_type.eq.3) then
6024 if(c2fit_unlike_1d(i) .ne. 0.0) then
6025 chisq_unlike_1d = chisq_unlike_1d + ((c2fit_unlike_1d(i)
6026 1 - c2mod_unlike_1d(i))/c2err_unlike_1d(i))**2
6030 end if ! End 1D correlation function, chi-square option
6032 if(switch_3d .gt. 0) then
6033 if(switch_type.eq.1 .or. switch_type.eq.3) then
6035 if(n_3d_fine .gt. 0) then
6039 if(c2fit_like_3d_fine(i,j,k).ne.0.0) then
6040 chisq_like_3d_fine = chisq_like_3d_fine
6041 1 + ((c2fit_like_3d_fine(i,j,k)
6042 2 - c2mod_like_3d_fine(i,j,k))
6043 3 /c2err_like_3d_fine(i,j,k))**2
6050 if(n_3d_coarse .gt. 0) then
6051 do i = 1,n_3d_coarse
6052 do j = 1,n_3d_coarse
6053 do k = 1,n_3d_coarse
6054 if((i+j+k).gt.3) then
6055 if(c2fit_like_3d_coarse(i,j,k).ne.0.0) then
6056 chisq_like_3d_coarse = chisq_like_3d_coarse
6057 1 +((c2fit_like_3d_coarse(i,j,k)
6058 2 - c2mod_like_3d_coarse(i,j,k))
6059 3 /c2err_like_3d_coarse(i,j,k))**2
6069 if(switch_type.eq.2 .or. switch_type.eq.3) then
6071 if(n_3d_fine .gt. 0) then
6075 if(c2fit_unlike_3d_fine(i,j,k).ne.0.0) then
6076 chisq_unlike_3d_fine = chisq_unlike_3d_fine
6077 1 + ((c2fit_unlike_3d_fine(i,j,k)
6078 2 - c2mod_unlike_3d_fine(i,j,k))
6079 3 /c2err_unlike_3d_fine(i,j,k))**2
6086 if(n_3d_coarse .gt. 0) then
6087 do i = 1,n_3d_coarse
6088 do j = 1,n_3d_coarse
6089 do k = 1,n_3d_coarse
6090 if((i+j+k).gt.3) then
6091 if(c2fit_unlike_3d_coarse(i,j,k).ne.0.0) then
6092 chisq_unlike_3d_coarse = chisq_unlike_3d_coarse
6093 1 +((c2fit_unlike_3d_coarse(i,j,k)
6094 2 - c2mod_unlike_3d_coarse(i,j,k))
6095 3 /c2err_unlike_3d_coarse(i,j,k))**2
6104 end if ! End of 3D Correlation Function, Chi-Square Option
6106 CCC Obtain chi-squares for one-body distributions
6108 if(pid(1) .gt. 0) then
6109 n1fac = float(n_part_used_1_trk)/float(n_part_used_1_ref)
6110 avgerrsq_pt_1 = float(n_part_used_1_trk)/float(n_pt_bins)
6111 avgerrsq_phi_1 = float(n_part_used_1_trk)/float(n_phi_bins)
6112 avgerrsq_eta_1 = float(n_part_used_1_trk)/float(n_eta_bins)
6113 avgerrsq_pt_ref_1 = float(n_part_used_1_ref)/float(n_pt_bins)
6114 avgerrsq_phi_ref_1 = float(n_part_used_1_ref)/float(n_phi_bins)
6115 avgerrsq_eta_ref_1 = float(n_part_used_1_ref)/float(n_eta_bins)
6118 if(pid(2) .gt. 0) then
6119 n2fac = float(n_part_used_2_trk)/float(n_part_used_2_ref)
6120 avgerrsq_pt_2 = float(n_part_used_2_trk)/float(n_pt_bins)
6121 avgerrsq_phi_2 = float(n_part_used_2_trk)/float(n_phi_bins)
6122 avgerrsq_eta_2 = float(n_part_used_2_trk)/float(n_eta_bins)
6123 avgerrsq_pt_ref_2 = float(n_part_used_2_ref)/float(n_pt_bins)
6124 avgerrsq_phi_ref_2 = float(n_part_used_2_ref)/float(n_phi_bins)
6125 avgerrsq_eta_ref_2 = float(n_part_used_2_ref)/float(n_eta_bins)
6128 if(pid(1) .gt. 0) then
6129 if(mode .eq. 1) then
6132 1 chisq1(hist1_pt_1,href1_pt_1,max_h_1d,avgerrsq_pt_1,
6133 2 avgerrsq_pt_ref_1,n1fac,n_pt_bins)
6134 3 +chisq1(hist1_phi_1,href1_phi_1,max_h_1d,avgerrsq_phi_1,
6135 4 avgerrsq_phi_ref_1,n1fac,n_phi_bins)
6136 5 +chisq1(hist1_eta_1,href1_eta_1,max_h_1d,avgerrsq_eta_1,
6137 6 avgerrsq_eta_ref_1,n1fac,n_eta_bins)
6139 else if(mode .eq. 2) then
6142 1 chisq1(htmp1_pt_1,href1_pt_1,max_h_1d,avgerrsq_pt_1,
6143 2 avgerrsq_pt_ref_1,n1fac,n_pt_bins)
6144 3 +chisq1(htmp1_phi_1,href1_phi_1,max_h_1d,avgerrsq_phi_1,
6145 4 avgerrsq_phi_ref_1,n1fac,n_phi_bins)
6146 5 +chisq1(htmp1_eta_1,href1_eta_1,max_h_1d,avgerrsq_eta_1,
6147 6 avgerrsq_eta_ref_1,n1fac,n_eta_bins)
6150 end if ! End pid(1) one-body histogram chi-square calculation
6152 if(pid(2) .gt. 0) then
6153 if(mode .eq. 1) then
6156 1 chisq1(hist1_pt_2,href1_pt_2,max_h_1d,avgerrsq_pt_2,
6157 2 avgerrsq_pt_ref_2,n2fac,n_pt_bins)
6158 3 +chisq1(hist1_phi_2,href1_phi_2,max_h_1d,avgerrsq_phi_2,
6159 4 avgerrsq_phi_ref_2,n2fac,n_phi_bins)
6160 5 +chisq1(hist1_eta_2,href1_eta_2,max_h_1d,avgerrsq_eta_2,
6161 6 avgerrsq_eta_ref_2,n2fac,n_eta_bins)
6163 else if(mode .eq. 2) then
6166 1 chisq1(htmp1_pt_2,href1_pt_2,max_h_1d,avgerrsq_pt_2,
6167 2 avgerrsq_pt_ref_2,n2fac,n_pt_bins)
6168 3 +chisq1(htmp1_phi_2,href1_phi_2,max_h_1d,avgerrsq_phi_2,
6169 4 avgerrsq_phi_ref_2,n2fac,n_phi_bins)
6170 5 +chisq1(htmp1_eta_2,href1_eta_2,max_h_1d,avgerrsq_eta_2,
6171 6 avgerrsq_eta_ref_2,n2fac,n_eta_bins)
6174 end if ! End pid(2) one-body histogram chi-square calculation
6179 C----------------------------------------------------------------------
6182 real*4 function chisq1(h,href,maxh,herravgsq,hreferravgsq,
6186 CCC Compute chi-square for 1D histogram h(), with respect to the
6187 CCC reference histogram, href().
6189 C h(maxh) = 1D histogram array
6190 C href(maxh) = 1D reference histogram array
6191 C maxh = dimension of histogram arrays
6192 C herravgsq = average error squared in histogram h's bins
6193 C hreferravgsq = average error squared in ref. hist. href's bins
6194 C numfac = ratio of total number of entries in h to that
6196 C nbins = # bins to use in chi-square sum, starting at array
6197 C element 1,2,... nbins (where nbins .le. maxh)
6199 C The chi-square value is returned in chisq1
6201 CCC Local Variable Type Declarations:
6203 integer*4 maxh, nbins, i
6204 integer*4 h(maxh),href(maxh)
6206 real*4 herravgsq,hreferravgsq,numfac,numfacsq
6207 real*4 herrsq,hreferrsq
6210 numfacsq = numfac*numfac
6213 if(h(i) .gt. 0) then
6214 herrsq = float(h(i))
6219 if(href(i) .gt. 0) then
6220 hreferrsq = float(href(i))
6222 hreferrsq = hreferravgsq
6225 chisq1 = chisq1 + ((float(h(i)) - numfac*float(href(i)))**2)
6226 1 /(herrsq + numfacsq*hreferrsq)
6232 C-----------------------------------------------------------------------
6235 Subroutine write_data(mode,ievent)
6238 CCC This subroutine writes the main output file, 'hbt_simulation.out'
6239 C on File Unit 8. File Unit 8 is opened and closed by the main
6242 C Also, the computed 1- and 2-body reference histograms are printed
6243 C out from this subroutine on File Units 11 and 9, respectively. These
6244 C files are opened/closed here.
6246 C Output content determined by input parameter 'mode', where:
6248 C Mode Description of Output
6249 C ----- -----------------------------------------------------------
6250 C 1 basic output file header
6251 C input and derived quantities
6253 C 2 reference histograms (1 and 2-body)
6254 C saved to separate I/O File Unit=11,9 respectively
6256 C 3 reference histogram output
6258 C 4 correlation model
6260 C 5 correlation fit and one-body distributions
6261 C for each event, optional output
6263 C 6 inclusive one-body distributions and inclusive
6264 C correlation fit; projection onto 1D axes.
6267 Include 'common_parameters.inc'
6268 Include 'common_mesh.inc'
6269 Include 'common_histograms.inc'
6270 Include 'common_correlations.inc'
6271 Include 'common_coulomb.inc'
6272 Include 'common_event_summary.inc'
6274 Include 'common_track.inc'
6275 Include 'common_sec_track.inc'
6276 Include 'common_sec_track2.inc'
6278 CCC Local Variable Type Declarations:
6280 integer*4 mode,i,j,k,ievent,ref_print,nev
6282 real*4 nfac1,nfac2,ref_error
6283 real*4 c2mod_proj1(max_c2_3d)
6284 real*4 c2mod_proj2(max_c2_3d)
6285 real*4 c2mod_proj3(max_c2_3d)
6286 real*4 c2fit_proj1(max_c2_3d)
6287 real*4 c2fit_proj2(max_c2_3d)
6288 real*4 c2fit_proj3(max_c2_3d)
6289 real*4 c2err_proj1(max_c2_3d)
6290 real*4 c2err_proj2(max_c2_3d)
6291 real*4 c2err_proj3(max_c2_3d)
6293 C-------------------------------------------
6294 If(mode.eq.1) Then !Basic Output Header
6295 C-------------------------------------------
6300 C write(8,102) n_events
6301 write(8,103) n_pid_types,pid(1),mass1,pid(2),mass2
6302 write(8,104) ref_control
6303 write(8,105) switch_1d
6304 write(8,106) switch_3d
6305 write(8,107) switch_type
6306 write(8,108) switch_coherence
6307 write(8,109) switch_coulomb
6308 write(8,110) switch_fermi_bose
6309 write(8,1101) trk_accep
6310 write(8,111) print_full,print_sector_data
6311 C write(8,112) n_part_used_1_ref,n_part_used_2_ref
6312 C write(8,113) n_part_used_1_inc,n_part_used_2_inc
6313 C write(8,114) num_pairs_like_ref,num_pairs_unlike_ref
6314 C write(8,115) num_pairs_like_inc,num_pairs_unlike_inc
6317 write(8,118) Rside,Rout,Rlong
6318 write(8,119) Rperp,Rparallel,R0
6324 write(8,125) chisq_wt_like_1d
6325 write(8,126) chisq_wt_unlike_1d
6326 write(8,127) chisq_wt_like_3d_fine
6327 write(8,128) chisq_wt_unlike_3d_fine
6328 write(8,129) chisq_wt_like_3d_coarse
6329 write(8,130) chisq_wt_unlike_3d_coarse
6330 write(8,131) chisq_wt_hist1_1
6331 write(8,132) chisq_wt_hist1_2
6333 write(8,134) n_pt_bins,pt_bin_size,pt_min,pt_max
6334 write(8,135) n_phi_bins,phi_bin_size,phi_min,phi_max
6335 write(8,136) n_eta_bins,eta_bin_size,eta_min,eta_max
6337 write(8,138) n_px_bins,delpx,px_min,px_max
6338 write(8,139) n_py_bins,delpy,py_min,py_max
6339 write(8,140) n_pz_bins,delpz,pz_min,pz_max
6340 write(8,141) n_sectors
6342 write(8,143) n_1d_fine,n_1d_coarse,n_1d_total
6343 write(8,144) binsize_1d_fine,binsize_1d_coarse
6344 write(8,145) qmid_1d,qmax_1d
6346 write(8,147) n_3d_fine,n_3d_coarse,n_3d_total
6347 write(8,148) binsize_3d_fine,binsize_3d_coarse
6348 write(8,149) qmid_3d,qmax_3d
6349 write(8,150) n_3d_fine_project
6351 CCC Formats for Mode=1 Output
6353 100 format( 15x,50('*'))
6354 101 format( 15x,'*****',7x,'HBT CORRELATION SIMULATION',7x,'*****')
6355 102 format(///15x,'Number of Events in Event Text Input File=',I5)
6356 103 format( /15x,'#PID types=',I2,' PID#,mass=',I2,F8.5,
6357 1' PID#,mass=',I2,F8.5)
6358 104 format(/ 15x,'Reference Spectra Selection Option=',I2)
6359 105 format(// 15x,'Control Switches: Switch_1d =',I2)
6360 106 format( 15x,' Switch_3d =',I2)
6361 107 format( 15x,' Switch_type =',I2)
6362 108 format( 15x,' Switch_coherence =',I2)
6363 109 format( 15x,' Switch_coulomb =',I2)
6364 110 format( 15x,' Switch_fermi_bose =',I2)
6365 1101 format( 15x,' trk_accep =',F10.7)
6366 111 format(/ 15x,'Print Options: Full=',I2,' Sectors=',I2)
6368 1'Number particles used in Reference, for PID types=',2I5)
6370 1'Number particles used in Inclusive, for PID types=',2I5)
6372 1'Number pairs used in Reference, like and unlike=',2I5)
6374 1'Number pairs used in Inclusive, like and unlike=',2I5)
6375 116 format(// 15x,'Correlation Model Parameters: Chaoticity =',F8.5)
6376 117 format( 15x,'1D Spherical Source Radius=',F8.4)
6377 118 format( 15x,'Bertsch-Pratt R-side,out,long=',3F8.4)
6378 119 format( 15x,'YKP R-perp,parallel,time=',3F8.4)
6379 120 format( 15x,'Coulomb parameter=',F8.5)
6380 121 format(// 15x,'Iteration Controls: Random # seed =',I10)
6381 122 format( 15x,' Max # iterations =',I5)
6382 123 format( 15x,' Momentum Shift Range =',F8.5)
6383 124 format( 15x,' Min % Chi-Sq limit =',F8.5)
6384 125 format(// 15x,'CHI-Sq Weights: correl,like, 1d =',F8.5)
6385 126 format( 15x,' correl,unlike,1d =',F8.5)
6386 127 format( 15x,' correl,like, 3d_fine =',F8.5)
6387 128 format( 15x,' correl,unlike,3d_fine =',F8.5)
6388 129 format( 15x,' correl,like, 3d_coarse=',F8.5)
6389 130 format( 15x,' correl,unlike,3d_coarse=',F8.5)
6390 131 format( 15x,' 1-body, PID#1 =',F8.5)
6391 132 format( 15x,' 1-body, PID#2 =',F8.5)
6392 133 format(// 15x,'Momentum Space Acceptance Range and 1D Bins:')
6393 134 format( 15x,'#bins,bin size,min,max for pt =',I5,3F8.5)
6394 135 format( 15x,'#bins,bin size,min,max for phi=',I5,3F8.4)
6395 136 format( 15x,'#bins,bin size,min,max for eta=',I5,3F8.4)
6396 137 format(// 15x,'Momentum Space Sectors:')
6397 138 format( 15x,'#sectors,sectorsize,min,max for px=',I5,3F8.4)
6398 139 format( 15x,'#sectors,sectorsize,min,max for py=',I5,3F8.4)
6399 140 format( 15x,'#sectors,sectorsize,min,max for pz=',I5,3F8.4)
6400 141 format( 15x,'Total Number of Sectors =',I5)
6401 142 format(// 15x,'2-Body Correlations, 1-D Grid:')
6402 143 format( 15x,'#bins_fine,coarse,total =',3I5)
6403 144 format( 15x,'bin size - fine, coarse =',2F8.5)
6404 145 format( 15x,'Q mid point, Q maximum =',2F8.5)
6405 146 format(// 15x,'2-Body Correlations, 3-D Grid:')
6406 147 format( 15x,'#bins - fine, coarse, total =',3I5)
6407 148 format( 15x,'bin size - fine, coarse =',2F8.5)
6408 149 format( 15x,'Q mid point, Q maximum =',2F8.5)
6409 150 format( 15x,'# 3D fine bin projected =',I5)
6411 CCC END mode=1 Output and Formats
6413 C-----------------------------
6414 Else If(mode.eq.2) Then !Store 2- and 1-body Ref. Histograms
6415 C-----------------------------
6418 open(unit=9,status='unknown',access='sequential',
6419 1 name='hbt_pair_reference.hist')
6420 open(unit=11,status='unknown',access='sequential',
6421 1 name='hbt_singles_reference.hist')
6423 C Write Pair Reference Hist:
6425 write(9,201) n_pid_types,pid(1),pid(2)
6426 write(9,202) n_pt_bins,pt_min,pt_max
6427 write(9,202) n_phi_bins,phi_min,phi_max
6428 write(9,202) n_eta_bins,eta_min,eta_max
6429 write(9,201) switch_1d,switch_3d,switch_type
6430 write(9,203) n_1d_fine,n_1d_coarse,n_3d_fine,n_3d_coarse
6431 write(9,204) binsize_1d_fine,binsize_1d_coarse,
6432 1 binsize_3d_fine,binsize_3d_coarse
6433 write(9,201) num_pairs_like_ref,num_pairs_unlike_ref
6435 202 format(2x,I10,2E15.6)
6437 204 format(2x,4E15.6)
6440 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6441 if(switch_type.eq.1.or.switch_type.eq.3) then
6442 write(9,205) (href_like_1d(i),i=1,n_1d_total)
6444 if(switch_type.eq.2.or.switch_type.eq.3) then
6445 write(9,205) (href_unlike_1d(i),i=1,n_1d_total)
6447 endif !End 1D Ref. Hist. Output
6449 if(switch_3d.gt.0) then
6450 if(switch_type.eq.1.or.switch_type.eq.3) then
6452 if(n_3d_fine.gt.0) then
6456 write(9,205) href_like_3d_fine(i,j,k)
6462 if(n_3d_coarse.gt.0) then
6466 write(9,205) href_like_3d_coarse(i,j,k)
6474 if(switch_type.eq.2.or.switch_type.eq.3) then
6476 if(n_3d_fine.gt.0) then
6480 write(9,205) href_unlike_3d_fine(i,j,k)
6486 if(n_3d_coarse.gt.0) then
6490 write(9,205) href_unlike_3d_coarse(i,j,k)
6497 endif !End 3D Reference Histograms Output
6499 CC Write One-Body - singles histograms:
6501 write(11,201) n_pid_types,pid(1),pid(2)
6502 write(11,202) n_pt_bins,pt_min,pt_max
6503 write(11,202) n_phi_bins,phi_min,phi_max
6504 write(11,202) n_eta_bins,eta_min,eta_max
6505 write(11,201) n_part_used_1_ref,n_part_used_2_ref
6507 if(pid(1).gt.0) then
6508 write(11,205)(href1_pt_1(i),i=1,n_pt_bins)
6509 write(11,205)(href1_phi_1(i),i=1,n_phi_bins)
6510 write(11,205)(href1_eta_1(i),i=1,n_eta_bins)
6514 if(pid(2).gt.0) then
6515 write(11,205)(href1_pt_2(i),i=1,n_pt_bins)
6516 write(11,205)(href1_phi_2(i),i=1,n_phi_bins)
6517 write(11,205)(href1_eta_2(i),i=1,n_eta_bins)
6523 CCC END mode=2 Reference Histogram Output
6525 C-----------------------------
6526 Else If(mode.eq.3) Then !Print out the Reference Histograms
6527 C-----------------------------
6531 write(8,302) n_pt_bins,pt_min,pt_max
6532 write(8,303) n_phi_bins,phi_min,phi_max
6533 write(8,304) n_eta_bins,eta_min,eta_max
6534 write(8,305) n_part_used_1_ref,n_part_used_2_ref
6538 write(8,307) i,href1_pt_1(i),href1_pt_2(i)
6543 write(8,307) i,href1_phi_1(i),href1_phi_2(i)
6548 write(8,307) i,href1_eta_1(i),href1_eta_2(i)
6552 write(8,311) n_1d_fine,n_1d_coarse
6553 write(8,312) binsize_1d_fine,binsize_1d_coarse
6554 write(8,313) n_3d_fine,n_3d_coarse
6555 write(8,314) binsize_3d_fine,binsize_3d_coarse
6556 write(8,315) num_pairs_like_ref,num_pairs_unlike_ref
6558 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6561 write(8,307) i,href_like_1d(i),href_unlike_1d(i)
6563 endif !End Print Out of 2-body, 1D reference histogram
6565 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6570 write(8,318) i,j,k,href_like_3d_fine(i,j,k),
6571 1 href_unlike_3d_fine(i,j,k)
6575 endif !End Print Out of 2-Body, 3D-Fine Mesh Ref. Hist.
6578 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
6583 write(8,318) i,j,k,href_like_3d_coarse(i,j,k),
6584 1 href_unlike_3d_coarse(i,j,k)
6588 endif !End Print Out of 2-Body, 3D-Coarse Mesh Ref. Hist.
6590 CCC Formats for mode=3 Output:
6592 300 format(///5x,15('*'),'REFERENCE HISTOGRAMS',15('*'))
6593 301 format(//15x,'ONE-BODY REFERENCE DISTRIBUTIONS:')
6594 302 format(/ 15x,'PT BINS: (#,min,max)=',I5,2F8.4)
6595 303 format( 15x,'PHI BINS:(#,min,max)=',I5,2F8.4)
6596 304 format( 15x,'ETA BINS:(#,min,max)=',I5,2F8.4)
6597 305 format( 15x,'Number particles used in Ref, PID type1,2=',2I8)
6598 306 format(/ 9x,'PT',10x,'BIN#',5x,'PID-1',5x,'PID-2')
6599 308 format(/ 9x,'PHI',9x,'BIN#',5x,'PID-1',5x,'PID-2')
6600 309 format(/ 9x,'ETA',9x,'BIN#',5x,'PID-1',5x,'PID-2')
6601 307 format( 20x,I5,2I10)
6602 310 format(///15x,'TWO-BODY REFERENCE DISTRIBUTIONS:')
6603 311 format(/ 15x,'#BINS FOR 1D-Fine and Coarse Grid =',2I5)
6604 312 format( 15x,'BIN SIZES FOR 1D-Fine and Coarse =',2F8.5)
6605 313 format( 15x,'#BINS FOR 3D-Fine and Coarse Grid =',2I5)
6606 314 format( 15x,'BIN SIZES FOR 3D-Fine and Coarse =',2F8.5)
6607 315 format( 15x,'Number of Like and Unlike Pairs For Ref. = ',
6609 316 format(/5x,'2-BODY, 1D',6x,'BIN#',5x,'LIKE',5x,'UNLIKE')
6610 317 format(/3x,'2-BODY, 3D-FINE',2x,'BIN:i',4x,'j',4x,
6611 1 'k',5x,'LIKE',5x,'UNLIKE')
6612 318 format( 20x,3I5,2I10)
6613 319 format(/2x,'2-BODY, 3D-COARSE',1x,'BIN:i',4x,'j',4x,
6614 1 'k',5x,'LIKE',5x,'UNLIKE')
6616 CCC END mode=3 Output and Formats
6618 C----------------------------
6619 Else If(mode.eq.4) Then !Print Correlation Function Model
6620 C----------------------------
6623 write(8,311) n_1d_fine,n_1d_coarse
6624 write(8,312) binsize_1d_fine,binsize_1d_coarse
6625 write(8,313) n_3d_fine,n_3d_coarse
6626 write(8,314) binsize_3d_fine,binsize_3d_coarse
6629 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6632 write(8,407) i,c2mod_like_1d(i),c2mod_unlike_1d(i)
6634 endif !End Print Out of 2-body, 1D Model Correction Functions
6636 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6641 write(8,418) i,j,k,c2mod_like_3d_fine(i,j,k),
6642 1 c2mod_unlike_3d_fine(i,j,k)
6646 endif !End Print Out of 2-Body, 3D-Fine mesh Model Correl. Function
6648 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
6653 write(8,418) i,j,k,c2mod_like_3d_coarse(i,j,k),
6654 1 c2mod_unlike_3d_coarse(i,j,k)
6658 endif !End Print Out of 2-Body, 3D-Coarse Model Correlation Function
6660 if(switch_coulomb.eq.3) then ! Print interpolated Pratt Model
6661 CC ! Coulomb Correction for Finite
6662 CC ! Source Radius Q0
6666 write(8,403) i,q_coul(i),c2_coul_like(i),
6671 CCC Additional Formats for Mode=4 Output:
6673 400 format(///5x,15('*'),'MODEL CORRELATION FUNCTIONS',15('*'))
6674 401 format(///15x,'COULOMB SOURCE RADIUS FOR PRATT MODEL=',F8.4)
6675 402 format(// 5x,'q-bin',2x,'q',4x,'C2_coul_like',2x,
6677 403 format( 5x,I5,3E15.6)
6678 407 format( 20x,I5,2F10.7)
6679 418 format( 20x,3I5,2F10.7)
6681 CCC END MODE = 4 OUTPUT AND FORMATS
6683 C------------------------------
6684 Else If(mode.eq.5) Then ! Optional Output for 1- and 2-Body Fits
6685 C------------------------------ ! for each event.
6686 C write(*,*)'Event ', ievent
6687 C write(*,*)'chisq_total_store(event) ',chisq_total_store(ievent)
6689 write(8,501) n_part_1_trk,n_part_2_trk,n_part_tot_trk
6690 write(8,502) n_part_used_1_trk, n_part_used_2_trk
6691 write(8,503) num_pairs_like, num_pairs_unlike
6693 CCC Output one-body distributions for event:
6696 if(pid(1) .gt. 0) then
6697 nfac1 = float(n_part_used_1_trk)/float(n_part_used_1_ref)
6702 ref_print = int(nfac1*float(href1_pt_1(i)))
6703 ref_error = nfac1*sqrt(float(href1_pt_1(i)))
6704 write(8,510) i,hist1_pt_1(i),ref_print,ref_error
6709 ref_print = int(nfac1*float(href1_phi_1(i)))
6710 ref_error = nfac1*sqrt(float(href1_phi_1(i)))
6711 write(8,510) i,hist1_phi_1(i),ref_print,ref_error
6716 ref_print = int(nfac1*float(href1_eta_1(i)))
6717 ref_error = nfac1*sqrt(float(href1_eta_1(i)))
6718 write(8,510) i,hist1_eta_1(i),ref_print,ref_error
6721 end if ! End PID # 1, One-Body Distribution Output
6723 if(pid(2) .gt. 0) then
6724 nfac2 = float(n_part_used_2_trk)/float(n_part_used_2_ref)
6729 ref_print = int(nfac2*float(href1_pt_2(i)))
6730 ref_error = nfac2*sqrt(float(href1_pt_2(i)))
6731 write(8,510) i,hist1_pt_2(i),ref_print,ref_error
6736 ref_print = int(nfac2*float(href1_phi_2(i)))
6737 ref_error = nfac2*sqrt(float(href1_phi_2(i)))
6738 write(8,510) i,hist1_phi_2(i),ref_print,ref_error
6743 ref_print = int(nfac2*float(href1_eta_2(i)))
6744 ref_error = nfac2*sqrt(float(href1_eta_2(i)))
6745 write(8,510) i,hist1_eta_2(i),ref_print,ref_error
6748 end if ! End PID # 2, One-Body Distribution Output
6750 CCC Output Two-Body Correlation Functions for Event:
6753 if(switch_1d.gt.0 .and. n_1d_total.gt.0) then
6758 write(8,523) i,c2mod_like_1d(i),c2fit_like_1d(i),
6759 1 c2err_like_1d(i),c2mod_unlike_1d(i),
6760 2 c2fit_unlike_1d(i),c2err_unlike_1d(i)
6762 end if ! End 1D Correlation Model and Fit Output
6764 if(switch_3d.gt.0 .and. n_3d_fine.gt.0) then
6771 write(8,526) i,j,k,c2mod_like_3d_fine(i,j,k),
6772 1 c2fit_like_3d_fine(i,j,k),c2err_like_3d_fine(i,j,k),
6773 2 c2mod_unlike_3d_fine(i,j,k),c2fit_unlike_3d_fine(i,j,k),
6774 3 c2err_unlike_3d_fine(i,j,k)
6778 end if ! End 3D Fine Mesh Correlation Model and Fit Output
6780 if(switch_3d.gt.0 .and. n_3d_coarse.gt.0) then
6784 do i = 1,n_3d_coarse
6785 do j = 1,n_3d_coarse
6786 do k = 1,n_3d_coarse
6787 write(8,526) i,j,k,c2mod_like_3d_coarse(i,j,k),
6788 1 c2fit_like_3d_coarse(i,j,k),c2err_like_3d_coarse(i,j,k),
6789 2 c2mod_unlike_3d_coarse(i,j,k),c2fit_unlike_3d_coarse(i,j,k),
6790 3 c2err_unlike_3d_coarse(i,j,k)
6794 end if ! End 3D Coarse Mesh Correlation Model and Fit Output
6796 CCC Output Event Summary and Chi-Square Information for Event:
6799 write(8,540) num_iter(ievent)
6800 write(8,541) n_part_used_1_store(ievent),
6801 1 n_part_used_2_store(ievent)
6802 write(8,5411) n_part_tot_store(ievent)
6803 write(8,542) num_sec_flagged_store(ievent)
6804 write(8,543) frac_trks_out(ievent),frac_trks_flag(ievent)
6805 write(8,544) chisq_like_1d_store(ievent),
6806 1 chisq_unlike_1d_store(ievent)
6807 write(8,545) chisq_like_3d_fine_store(ievent),
6808 1 chisq_unlike_3d_fine_store(ievent)
6809 write(8,546) chisq_like_3d_coarse_store(ievent),
6810 1 chisq_unlike_3d_coarse_store(ievent)
6811 write(8,547) chisq_hist1_1_store(ievent),
6812 1 chisq_hist1_2_store(ievent)
6813 write(8,548) chisq_total_store(ievent)
6815 CCC Formats for Mode = 5 Output:
6817 500 Format(///5x,5('*'),'Fitted 1-Body Distributions and ',
6818 1 'Correlations for Event #',I5,5('*'))
6819 501 Format(//15x,'Number of Particles of PID types 1,2,total = ',
6821 502 Format( 15x,'Number of Particles of PID types 1,2 Used = ',
6823 503 Format( 15x,'Number of Pairs Used - Like and Unlike = ',2I10)
6824 504 Format(//5x,'Fitted and Normalized Reference One-Body ',
6826 505 Format( /10x,'Particle Type 1: Reference Scale Factor = ',E12.5)
6827 506 Format( /10x,'Particle Type 2: Reference Scale Factor = ',E12.5)
6828 507 Format(/2x,' PT: BIN#',5x,'hist1',7x,'href1-scaled',2x,
6830 508 Format(/2x,'PHI: BIN#',5x,'hist1',7x,'href1-scaled',2x,
6832 509 Format(/2x,'ETA: BIN#',5x,'hist1',7x,'href1-scaled',2x,
6834 510 Format(7x,I4,3x,I7,8x,I7,7x,F10.5)
6835 520 Format(//5x,'Model and Fitted Correlations')
6836 530 Format(//21x,'One-Dimensional Fit - Fine & Coarse Mesh')
6837 531 Format(//25x,'Three-Dimensional Fit - Fine Mesh')
6838 532 Format(//24x,'Three-Dimensional Fit - Coarse Mesh')
6839 521 Format(/1x,'BIN',13x,'LIKE PAIRS',27x,'UNLIKE PAIRS')
6840 522 Format(8x,'MOD',9x,'FIT',9x,'ERR',11x,'MOD',9x,'FIT',9x,
6842 523 Format(1x,I3,3E12.4,2x,3E12.4)
6843 524 Format(/2x,'BINS',12x,'LIKE PAIRS',24x,'UNLIKE PAIRS')
6844 525 Format(1x,' i j k',4x,'MOD',8x,'FIT',8x,'ERR',10x,'MOD',
6845 1 8x,'FIT',8x,'ERR',/)
6846 526 Format(1x,3I2,3E11.4,2x,3E11.4)
6847 539 Format(///10x,'Event and Chi-Square Summary for Event #',I5)
6848 540 Format( //15x,'Number of Iterations =',F10.2)
6849 541 Format( 15x,'# Particles Used for PID Types1,2=',2F10.2)
6850 5411 Format( 15x,'Total # Particles in track table =',F10.2)
6851 542 Format( 15x,'# Sectors Flagged =',F10.2)
6852 543 Format( 15x,'Frac Trks Out of Accep., Flagged =',2E11.4)
6853 544 Format( 15x,'Chi-Sq: 1D - Like & Unlike =',2E11.4)
6854 545 Format( 15x,'Chi-Sq: 3D - Fine -Like & Unlike =',2E11.4)
6855 546 Format( 15x,'Chi-Sq: 3D - Coarse-Like &Unlike =',2E11.4)
6856 547 Format( 15x,'Chi-Sq: One-Body Dist. PID# 1&2 =',2E11.4)
6857 548 Format( 15x,'Chi-Sq: Total Weighted =',E11.4)
6859 CCC End Mode = 5 Output and Formats
6861 C------------------------------
6862 Else If(mode.eq.6) Then ! Inclusive 1 & 2 Body Output
6863 C------------------------------
6864 write(8,600) n_events
6865 write(8,601) n_part_used_1_inc,n_part_used_2_inc
6866 write(8,602) num_pairs_like_inc,num_pairs_unlike_inc
6869 if(pid(1).gt.0) then
6870 C Division by zero check
6871 IF (n_part_used_1_ref .LE. 0) THEN
6872 PRINT*,'************************************'
6873 PRINT*,'* HBT PROCESSOR *'
6874 PRINT*,'* Number of particles selected for *'
6875 PRINT*,'* processing is less or equal *'
6876 PRINT*,'* !!!!!!!! ZER0 !!!!!!!!!! *'
6877 PRINT*,'* unable to proceed *'
6878 PRINT*,'* EXITING FORTRAN *'
6880 PRINT*,'* HINT: broad the parameter regions*'
6881 PRINT*,'* OR/AND number of particles OR/AND*'
6882 PRINT*,'* number of events *'
6883 PRINT*,'************************************'
6885 5481 FORMAT(5x,'Number of particles selected for processing is',
6886 1 ' less or equal 0',
6891 nfac1=float(n_part_used_1_inc)/float(n_part_used_1_ref)
6895 ref_print=int(nfac1*float(href1_pt_1(i)))
6896 ref_error=nfac1*sqrt(float(href1_pt_1(i)))
6897 write(8,510) i,hinc1_pt_1(i),ref_print,ref_error
6902 ref_print=int(nfac1*float(href1_phi_1(i)))
6903 ref_error=nfac1*sqrt(float(href1_phi_1(i)))
6904 write(8,510) i,hinc1_phi_1(i),ref_print,ref_error
6909 ref_print=int(nfac1*float(href1_eta_1(i)))
6910 ref_error=nfac1*sqrt(float(href1_eta_1(i)))
6911 write(8,510) i,hinc1_eta_1(i),ref_print,ref_error
6914 endif !END PID #1 One-BODY INCL. DISTRIBUTION OUTPUT
6916 if(pid(2).gt.0) then
6917 nfac2=float(n_part_used_2_inc)/float(n_part_used_2_ref)
6921 ref_print=int(nfac2*float(href1_pt_2(i)))
6922 ref_error=nfac2*sqrt(float(href1_pt_2(i)))
6923 write(8,510) i,hinc1_pt_2(i),ref_print,ref_error
6928 ref_print=int(nfac2*float(href1_phi_2(i)))
6929 ref_error=nfac2*sqrt(float(href1_phi_2(i)))
6930 write(8,510) i,hinc1_phi_2(i),ref_print,ref_error
6935 ref_print=int(nfac2*float(href1_eta_2(i)))
6936 ref_error=nfac2*sqrt(float(href1_eta_2(i)))
6937 write(8,510) i,hinc1_eta_2(i),ref_print,ref_error
6940 endif !END PID #2 One-BODY INCl. DISTRIBUTION OUTPUT
6942 CC OUTPUT TWO-BODY INCLUSIVE HISTOGRAMS:
6945 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6948 write(8,307) i,hinc_like_1d(i),hinc_unlike_1d(i)
6950 end if ! End Print out of 2-Body, 1D Inclusive Histograms
6952 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6957 write(8,318) i,j,k,hinc_like_3d_fine(i,j,k),
6958 1 hinc_unlike_3d_fine(i,j,k)
6962 endif ! End Print out of 2-Body, 3D-Fine Inclusive Histograms
6964 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
6969 write(8,318) i,j,k,hinc_like_3d_coarse(i,j,k),
6970 1 hinc_unlike_3d_coarse(i,j,k)
6974 endif ! End Print out of 2-Body, 3D-Coarse Inclusive Histograms
6976 CC OUTPUT TWO-BODY INCL.CORRELATION FUNCTIONS FOR EVENT
6979 if(switch_1d.gt.0.and.n_1d_total.gt.0) then
6984 write(8,523) i,c2mod_like_1d(i),c2fit_like_1d(i),
6985 1 c2err_like_1d(i),c2mod_unlike_1d(i),
6986 2 c2fit_unlike_1d(i),c2err_unlike_1d(i)
6990 if(switch_3d.gt.0.and.n_3d_fine.gt.0) then
6997 write(8,526) i,j,k,c2mod_like_3d_fine(i,j,k),
6998 1 c2fit_like_3d_fine(i,j,k),
6999 2 c2err_like_3d_fine(i,j,k),
7000 3 c2mod_unlike_3d_fine(i,j,k),
7001 4 c2fit_unlike_3d_fine(i,j,k),
7002 5 c2err_unlike_3d_fine(i,j,k)
7009 if(switch_3d.gt.0.and.n_3d_coarse.gt.0) then
7016 write(8,526) i,j,k,c2mod_like_3d_coarse(i,j,k),
7017 1 c2fit_like_3d_coarse(i,j,k),
7018 2 c2err_like_3d_coarse(i,j,k),
7019 3 c2mod_unlike_3d_coarse(i,j,k),
7020 4 c2fit_unlike_3d_coarse(i,j,k),
7021 5 c2err_unlike_3d_coarse(i,j,k)
7028 CCC Compute and Print 1D projections of 3D fine mesh C2 model,
7029 CCC fit and errors for like and unlike pairs.
7031 if(switch_3d .gt. 0 .and. n_3d_fine .gt. 0) then
7032 if(switch_type .eq. 1 .or. switch_type .eq. 3) then
7033 Call c2_3d_projected(hinc_like_3d_fine,
7034 1 href_like_3d_fine,c2mod_like_3d_fine,
7035 2 c2mod_proj1,c2mod_proj2,c2mod_proj3,
7036 3 c2fit_proj1,c2fit_proj2,c2fit_proj3,
7037 4 c2err_proj1,c2err_proj2,c2err_proj3,
7038 5 max_h_3d, max_c2_3d, n_3d_fine,
7039 6 n_3d_fine_project,num_pairs_like_inc,
7040 7 num_pairs_like_ref)
7045 write(8,658) i,c2mod_proj1(i),c2fit_proj1(i),c2err_proj1(i)
7050 write(8,658) i,c2mod_proj2(i),c2fit_proj2(i),c2err_proj2(i)
7055 write(8,658) i,c2mod_proj3(i),c2fit_proj3(i),c2err_proj3(i)
7057 end if ! End Like pair output
7059 if(switch_type .eq. 2 .or. switch_type .eq. 3) then
7060 Call c2_3d_projected(hinc_unlike_3d_fine,
7061 1 href_unlike_3d_fine,c2mod_unlike_3d_fine,
7062 2 c2mod_proj1,c2mod_proj2,c2mod_proj3,
7063 3 c2fit_proj1,c2fit_proj2,c2fit_proj3,
7064 4 c2err_proj1,c2err_proj2,c2err_proj3,
7065 5 max_h_3d, max_c2_3d, n_3d_fine,
7066 6 n_3d_fine_project,num_pairs_unlike_inc,
7067 7 num_pairs_unlike_ref)
7071 write(8,658) i,c2mod_proj1(i),c2fit_proj1(i),c2err_proj1(i)
7076 write(8,658) i,c2mod_proj2(i),c2fit_proj2(i),c2err_proj2(i)
7081 write(8,658) i,c2mod_proj3(i),c2fit_proj3(i),c2err_proj3(i)
7083 end if ! End Unlike pair output
7084 end if ! End 1D projections
7087 CCC EVENT AND CHISQ SUMMARY INFORMATION:
7089 if(n_events.le.max_events) then
7099 write(8,623) i,num_iter(i),n_part_used_1_store(i),
7100 1 n_part_used_2_store(i),
7101 2 num_sec_flagged_store(i),
7102 3 frac_trks_out(i),frac_trks_flag(i),
7103 4 chisq_total_store(i)
7106 write(8,6231) trk_maxlen
7109 write(8,6233) i,n_part_tot_store(i)
7114 write(8,625) i,chisq_like_1d_store(i),
7115 1 chisq_unlike_1d_store(i)
7121 write(8,625) i,chisq_like_3d_fine_store(i),
7122 1 chisq_unlike_3d_fine_store(i)
7128 write(8,625) i,chisq_like_3d_coarse_store(i),
7129 1 chisq_unlike_3d_coarse_store(i)
7135 write(8,625) i,chisq_hist1_1_store(i),
7136 1 chisq_hist1_2_store(i)
7139 CCC Output the Mean and RMS values for the Event Loop:
7143 write(8,631) niter_mean,niter_rms
7144 write(8,632) npart1_mean,npart1_rms
7145 write(8,633) npart2_mean,npart2_rms
7146 write(8,6331) npart_tot_mean, npart_tot_rms
7147 write(8,634) nsec_flag_mean,nsec_flag_rms
7148 write(8,635) frac_trks_out_mean,frac_trks_out_rms
7149 write(8,636) frac_trks_flag_mean,frac_trks_flag_rms
7150 write(8,637) chi_l1d_mean,chi_l1d_rms
7151 write(8,638) chi_u1d_mean,chi_u1d_rms
7152 write(8,639) chi_l3f_mean,chi_l3f_rms
7153 write(8,640) chi_u3f_mean,chi_u3f_rms
7154 write(8,641) chi_l3c_mean,chi_l3c_rms
7155 write(8,642) chi_u3c_mean,chi_u3c_rms
7156 write(8,643) chi_1_1_mean,chi_1_1_rms
7157 write(8,644) chi_1_2_mean,chi_1_2_rms
7158 write(8,645) chi_tot_mean,chi_tot_rms
7160 CCC FORMATS FOR MODE = 6 OUTPUT
7162 600 format(/// 2x,'FITTED 1-BODY DIST. AND CORRELATIONS ',
7163 1 'FOR INCLUSIVE SUM OF',I5,' EVENTS')
7164 601 format(// 15x,'Inclusive # Particles USED of PID ',
7166 602 format( 15x,'Inclusive # of pairs used; like/unlike=',
7168 603 format(// 5x,'Inclusive and Normalized Reference ',
7169 1 'One-Body Distributions')
7170 604 format(/ 10x,'Inclusive: Particle Type 1 - Reference ',
7171 1 'Scale Factor=',E12.5)
7172 605 format(/ 2x,' PT: BIN#',5x,'hinc1',7x,'href1-scaled',2x,
7174 606 format(/ 2x,'PHI: BIN#',5x,'hinc1',7x,'href1-scaled',2x,
7176 607 format(/ 2x,'ETA: BIN#',5x,'hinc1',7x,'href1-scaled',2x,
7178 608 format(/ 10x,'Inclusive: Particle Type 2 - ',
7179 1 'Reference Scale Factor=',E12.5)
7180 620 format(// 5x,'MODEL AND INCLUSIVE FITTED CORRELATIONS')
7181 621 format(// 15x,'Event and Chi-Square Summary Lists')
7182 622 format(/ 3x,'event',2x,'#iter',3x,'#PID1',4x,'#PID2',3x,
7183 1 '#sec-flg',3x,'frac-out',4x,'frac-flg',3x,'CHISQ-TOT')
7184 623 format(3x,I5,2x,F6.0,2(1x,F8.0),1x,F9.0,
7185 1 2(1x,F11.8),1x,E11.4)
7186 6231 format(/5x,'Max# tracks allowed in track table = ',I8)
7187 6232 format(/5x,'event',4x,'Tot# trks')
7188 6233 format(5x,I5,F12.2)
7189 624 format(/5x,'event',4x,'CHI_l1d',8x,'CHI_u1d')
7190 626 format(/5x,'event',4x,'CHI_l3f',8x,'CHI_u3f')
7191 627 format(/5x,'event',4x,'CHI_l3c',8x,'CHI_u3c')
7192 628 format(/5x,'event',4x,'CHI_1-1',8x,'CHI_1-2')
7193 625 format(5x,I5,2E15.6)
7194 629 format(// 10x,'Event and Chi-Square Summary - ',
7195 1 'Mean and RMS Values')
7196 630 format(/ 14x,'Quantity',15x,'Mean',11x,'RMS')
7197 631 format( 5x,'Number of Iterations ',2E15.6)
7198 632 format( 5x,'#PID Type 1 ',2E15.6)
7199 633 format( 5x,'#PID Type 2 ',2E15.6)
7200 6331 format( 5x,'Tot # Tracks in Table ',2E15.6)
7201 634 format( 5x,'#Sectors Flagged ',2E15.6)
7202 635 format( 5x,'Frac. Trks Out of Accept. ',2E15.6)
7203 636 format( 5x,'Frac. Trks Flagged ',2E15.6)
7204 637 format( 5x,'CHISQ like 1D ',2E15.6)
7205 638 format( 5x,'CHISQ unlike 1D ',2E15.6)
7206 639 format( 5x,'CHISQ like 3D Fine ',2E15.6)
7207 640 format( 5x,'CHISQ unlike 3D Fine ',2E15.6)
7208 641 format( 5x,'CHISQ like 3D Coarse ',2E15.6)
7209 642 format( 5x,'CHISQ unlike 3D Coarse ',2E15.6)
7210 643 format( 5x,'CHISQ 1 Body #1 ',2E15.6)
7211 644 format( 5x,'CHISQ 1 Body #2 ',2E15.6)
7212 645 format( 5x,'CHISQ Total ',2E15.6)
7213 650 format(//10x ,'Inclusive Three-Dimensional Projected Fits -',
7215 651 format( /25x ,'Like Pairs - Axis #1 ')
7216 652 format( /25x ,'Like Pairs - Axis #2 ')
7217 653 format( /25x ,'Like Pairs - Axis #3 ')
7218 654 format( /25x ,'Unlike Pairs - Axis #1 ')
7219 655 format( /25x ,'Unlike Pairs - Axis #2 ')
7220 656 format( /25x ,'Unlike Pairs - Axis #3 ')
7221 657 format( 2x,'BIN#',3x,'Model',8x,'Fit',8x,'Error')
7222 658 format( 3x,I3,3E12.4)
7223 660 format(// 5x,'INCLUSIVE TWO-BODY HISTOGRAMS')
7225 CCC END MODE = 6 OUTPUT AND FORMATS
7234 C-----------------------------------------------------------------------
7237 subroutine c2_3d_projected(h,href,c2mod,
7238 1 c2mod_proj1,c2mod_proj2,c2mod_proj3,
7239 2 c2fit_proj1,c2fit_proj2,c2fit_proj3,
7240 3 c2err_proj1,c2err_proj2,c2err_proj3,
7241 4 maxh,maxc2,n,n_proj,num_pairs_sig,
7246 CCC This Subroutine computes the projected two-body correlation
7247 CCC function for 3D distributions - fine mesh only; for both the
7248 CCC correlation model (weighted with the reference histogram) and
7249 CCC the inclusive correlation fit.
7251 CCC Description of Input Variables in the Argument List:
7253 CCC h(maxh,maxh,maxh) = 3D fine mesh inclusive signal histog.
7254 CCC href(maxh,maxh,maxh) = 3D fine mesh inclusive background hist.
7255 CCC c2mod(maxc2,maxc2,maxc2) = 3D fine mesh correlation model
7256 CCC maxh = Dimension of 3D fine mesh histogram arrays
7257 CCC maxc2 = Dimension of 3D fine mesh correlation function arrays
7258 CCC n = Number of bins to use
7259 CCC n_proj = Number of bins to integrate in (i,j) to project onto (k)
7260 CCC num_pairs_sig = # pairs used in signal histogram
7261 CCC num_pairs_bkg = # pairs used in background histogram
7263 CCC Description of Output quantities:
7265 CCC c2mod_proj1,2,3(maxc2) = Reference histogram weighted 1D projections
7266 CCC of C2 model function along {1,2,3} axes.
7267 CCC c2fit_proj1,2,3(maxc2) = Fitted 3D correlation function projected
7268 CCC onto {1,2,3} axes.
7269 CCC c2err_proj1,2,3(maxc2) = Error in fitted 3D correlation function
7270 CCC projected onto {1,2,3} axes.
7272 CCC Local Variable Type Declarations:
7274 integer*4 maxh,maxc2,n,n_proj,num_pairs_sig,num_pairs_bkg
7275 integer*4 h(maxh,maxh,maxh),href(maxh,maxh,maxh)
7278 real*4 c2mod(maxc2,maxc2,maxc2)
7279 real*4 c2mod_proj1(maxc2),c2mod_proj2(maxc2),c2mod_proj3(maxc2)
7280 real*4 c2fit_proj1(maxc2),c2fit_proj2(maxc2),c2fit_proj3(maxc2)
7281 real*4 c2err_proj1(maxc2),c2err_proj2(maxc2),c2err_proj3(maxc2)
7282 real*4 a,a_error,b,b_error
7283 real*4 sum1n,sum1d,sum2n,sum2d,sum3n,sum3d
7285 CCC Initialize arrays to zero:
7288 c2mod_proj1(i) = 0.0
7289 c2mod_proj2(i) = 0.0
7290 c2mod_proj3(i) = 0.0
7291 c2fit_proj1(i) = 0.0
7292 c2fit_proj2(i) = 0.0
7293 c2fit_proj3(i) = 0.0
7294 c2err_proj1(i) = 0.0
7295 c2err_proj2(i) = 0.0
7296 c2err_proj3(i) = 0.0
7299 CCC Project Reference spectra (histogram) weighted model correlation:
7310 sum1n = sum1n + c2mod(i,j,k)*float(href(i,j,k))
7311 sum1d = sum1d + float(href(i,j,k))
7312 sum2n = sum2n + c2mod(k,i,j)*float(href(k,i,j))
7313 sum2d = sum2d + float(href(k,i,j))
7314 sum3n = sum3n + c2mod(j,k,i)*float(href(j,k,i))
7315 sum3d = sum3d + float(href(j,k,i))
7318 if(sum1d .le. 0.0) then
7319 c2mod_proj1(i) = 0.0
7321 c2mod_proj1(i) = sum1n/sum1d
7323 if(sum2d .le. 0.0) then
7324 c2mod_proj2(i) = 0.0
7326 c2mod_proj2(i) = sum2n/sum2d
7328 if(sum3d .le. 0.0) then
7329 c2mod_proj3(i) = 0.0
7331 c2mod_proj3(i) = sum3n/sum3d
7335 CCC Calculate and Project the fitted correlation functions:
7346 sum1n = sum1n + float(h(i,j,k))
7347 sum1d = sum1d + float(href(i,j,k))
7348 sum2n = sum2n + float(h(k,i,j))
7349 sum2d = sum2d + float(href(k,i,j))
7350 sum3n = sum3n + float(h(j,k,i))
7351 sum3d = sum3d + float(href(j,k,i))
7354 if(sum1n .le. 0.0 .or. sum1d .le. 0.0) then
7355 c2fit_proj1(i) = 0.0
7356 c2err_proj1(i) = 1.0
7358 a = sum1n/float(num_pairs_sig)
7359 a_error = sqrt(sum1n)/float(num_pairs_sig)
7360 b = sum1d/float(num_pairs_bkg)
7361 b_error = sqrt(sum1d)/float(num_pairs_bkg)
7362 c2fit_proj1(i) = a/b
7363 c2err_proj1(i) = c2fit_proj1(i)*sqrt((a_error/a)**2
7366 if(sum2n .le. 0.0 .or. sum2d .le. 0.0) then
7367 c2fit_proj2(i) = 0.0
7368 c2err_proj2(i) = 1.0
7370 a = sum2n/float(num_pairs_sig)
7371 a_error = sqrt(sum2n)/float(num_pairs_sig)
7372 b = sum2d/float(num_pairs_bkg)
7373 b_error = sqrt(sum2d)/float(num_pairs_bkg)
7374 c2fit_proj2(i) = a/b
7375 c2err_proj2(i) = c2fit_proj2(i)*sqrt((a_error/a)**2
7378 if(sum3n .le. 0.0 .or. sum3d .le. 0.0) then
7379 c2fit_proj3(i) = 0.0
7380 c2err_proj3(i) = 1.0
7382 a = sum3n/float(num_pairs_sig)
7383 a_error = sqrt(sum3n)/float(num_pairs_sig)
7384 b = sum3d/float(num_pairs_bkg)
7385 b_error = sqrt(sum3d)/float(num_pairs_bkg)
7386 c2fit_proj3(i) = a/b
7387 c2err_proj3(i) = c2fit_proj3(i)*sqrt((a_error/a)**2
7395 C----------------------------------------------------------------------
7397 C>>>>>>>>>>>>>> Piotr, this needs to be replaced with Ali random
7398 C>>>>>>>>>>>>>>> number generator
7400 * real*4 function hbtpran(i)
7404 * Call ranhbtp(r,1,i)
7409 * Include 'ranlux2.f'
7410 C----------------------------------------------------------------------