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