Welcome to the Phobos Multiplicity Reconstruction User's Guide
Robin Verdier
Updated March 18, 1999
Multiplicity reconstruction is carried out by a procedural controller class, TPhMul, whose methods
Back to top
To reconstruct multiplicities, you must have access to a complete copy of the compiled Phat environment, either a public version or a personal version that you or somebody else has created. It may be possible to access only the executable image and shared libraries, but this is not recommended. To use Phat, you must have defined the environment variables ROOTSYS, PHATHOME, PATH, and LD_LIBRARY_PATH correctly, as described in detail in Appendix A and in the Multiplicity Status and Multiplicity FAQ documents.
We will describe here ways to run the multiplicity and vertex reconstruction alone. Phat macros that run the entire online system may be found in $PHATHOME/macros. The following optional definition is convenient for the examples below:
setenv mm $PHATHOME/phatmul/macros
If this command fails, you are probably in the wrong shell; see Appendix A.
To reconstruct multiplicities, after defining the DISPLAY environment variable, type
phat $mm/runmul.C
If the phat command is not accepted, you must set your PATH variable as previously noted.
To run in batch mode, that is, without an X-display, type phat -b $mm/runmul.C
Respond to the dialogue questions about help, data file, data format, # of events to process and to skip, whether to compare the reconstruction with MC information if present, whether to use the standard vertex finder, and setup and output scripts, called macros in Phat. Press return to accept the default, or normal choice, for any option. Note that the default .root files are those accessible on ppcln; you must locate data files at the site on which you're running; see the Multiplicity Status document or Data Tracker documents for help.
You only need a setup macro if you want to reset parameter values such as the energy limit cut, vary geometry, tune, or display characters, turn on diagnostic histogtrams, or increase debugging report levels. To do this, copy $mm/setmul.C to your working directory, make it writeable using chmod +w setmul.C, and edit it using the given examples as a guide; then specify it as the setup macro when prompted by runmul.C.
You only need an output macro if you want to write the input data and results of the reconstruction into a DST for later plotting or reanalyzing. The standard choice, $mm/writemul.C, writes about 1 MB / event of compressed information containing all TPhMapSets plus additional plots and histograms. To make your own output choice, copy $mm/writemul.C to your working directory, make it writeable using chmod +w setmul.C, and edit it using the given examples as a guide; then name it as the output macro.
Once the events have been processed, you can plot results. There is a Phat macro for producing a standard set; to use it, type
.x $mm/plotmul.C
and follow the prompts. The standard display shows plots of occupancy, dE/track, total tracks/ primaries, and dN/dh for charged primaries, reading clockwise from the top left; as indicated on the legends, the blue markers indicate generated tracks and hits, and read as zeroes when no MC hits and tracks are available. The default plot is for the display names '*+in', which shows the combined results from all sensors ['*'], integrated over f [+in'] and averaged over events; other choices are possible, such as '*+in-1' for the most recent single event.
You can also create your own plots, using standard objects of class TGraphErrors or TH2F (histograms). To do thisyou need to know the names of the created TPhMul processor and of the Sensor maps and Physics displays, and the indices of the calculated variables. Macro runmul.C creates a TPhMul object named mp on the stack, so that its methods are invoked by the "." rather than the "->" notation.
To exhibit the names of sensors active for multiplicity reconstruction, enter
mp.ListSensors()
and, for active displays,
mp.ListDisplays()
The calculated variables are listed and defined in the header file $mm/../inc/TPhMul.h, q.v. The user interfaces are still evolving. For sensors, the useful variable indices are 0 for deposited energy, 1 for the number of MC tracks, and 2 for the number of MC primary tracks incident on each pad. The most important display indices 13 and 14, for the reconstructed multiplicity and its statistical error, and 23 and 24, for the generated MC equivalents. Other variables contain phase space weights, occupancy, weighted global corrections, and other variables intermediate to the reconstruction.
To make a plot, you first extract a standard ROOT plot object, either a 2-d histogram (class TH2F) or a graph with errors (class TGraphErrors). For definiteness, suppose you have processed some number of events and wish to look at reconstructed dN/dh. You would enter
d1 = mp.GetGraph("d1","*+in",13,14,360)
This statement uses the ROOT interpreter's shorthand; if you put it in a compiled C++ program, you would have to use the full C++ syntax:
TGraphErrors *d1 = mp.GetGraph("d1","*+in",13,14,360);
Here the symbol d1 is a pointer to a TGraphErrors object created by mp, while the text string "d1" is the name by which the object will be known to ROOT. Both are arbitrary, EXCEPT that if first letter of the ROOT name is "s" or "S", mp looks in its list of sensor maps, whereas for any other initial letter it looks in its list of displays.
Using wildcards to create combined displays can lead to poorly weighted results; the correct approach is to create a single display that contains the desired combination of detectors. This can be done dynamically.
The normalization, 360, converts the calculated d2N/(dhdf) to an integral over f. For the double-differential displays, the normalization should in general be 360 divided by the number of summed f-bins.
Other arguments allow you to choose the number of x- and y- bins in the plot, and the numbers of edge bins to skip (for dealing with small statistics).
Once you have extracted an object, you can draw it using its Draw methods, say
d1->Draw("ape").
The 'a' draws axes; 'p' puts a marker at each point, and 'e' shows error bars. You can add 'l' to draw a straight line between each pair of points. You should first create a TCanvas to draw on, but if you don't, ROOT will do it for you, and you can then resize it.
d1->Draw() gives you nothing. TGraph methods can be used to make the plots usable; for example,
d1->SetMarkerStyle(24)
d1->SetMarkerColor(5)
d1->GetHistogram()->SetXTitle("Eta")
d1->GetHistogram()->SetYTitle("Reconstructed dN/dEta")
d1->GetHistogram()->SetTitleOffset(1.4,"Y")
TCanvas provides a method that creates PostScript file for printing,
c1->Print("dNdEtaRecon.ps")
If the extension is .eps, the file will be in Encapsulated PostScript. These commands usually executed by a script, using .x as above. The script may contain [almost] any C++ code, but, in addition, must begin with a curly brace { and end with } . At present, ROOT does not allow a TGraph's axes to be labelled until the graph has actually been drawn, although the labelling for a TH2F can be specified on at creation. A request to the ROOT development team to correct this has not been answered.
You can get a lego plot of d2N/(dhdf) accumulated for all events with
d2 = mp.GetHist("d2","*-*1",13,14,360)
d2->Draw("lego2")
will give you color-coded heights.
To see the raw energy deposition in all the octagon sensors, you could do
d3=mp.GetHist("s1","MO*",0)
There are also tools to make formatted dumps of these objects:
.L $mm/printhist.C
printhist(d3)
and, similarly, printgraph for TGraphs.
It is your responsibility to delete the created objects, for example,
d3->Delete()
If you do not do so, they will gradually use up the heap and you could end up running out of memory. This is the reason that the plot methods return pointers. You can also get a pointer to any named object, such as "s1", by using ROOT's GetObject methods.
You can also access the information in each TPhPadMap directly; for example, the statement
m1 = mp.GetDisplay("*+in");
creates a pointer-to-TPhPadMap called m1 to point to the integrated all-events display map, so that
m1->Print(13,14);
shows the contents. You can get the same information with less typing:
mp.GetDisplay("*+in")->Print(13,14);
There are standard scripts to do more elaborate displays, such as $mm/plotmul.C.
Back to top
Multiplicity reconstruction is based on a procedural class, TPhMul, which implements the reconstruction algorithms and provides access to the methods of classes that TPhMul invokes, and several object classes. The user's basic interface to TPhMul is through the method Process(), which in turn calls methods GetData, to bring real or MC data into the sensor maps, and Recon(), which carries out the actual reconstruction. Other TPhMul methods, GetSensMap, GetDisplay, and GetDiag, return pointers to the stored data and physics variables that are necessary for displaying data. The object classes are described in the following section.
Back to top
TPhMul::Recon:
This method reconstructs the most probable number of charged primary tracks incident on each active sensor pad, by normalizing the deposited energy to normal incidence assuming that all tracks came from the interaction vertex, limiting the result to a specified cutoff value, normally 500 keV, and dividing by the average energy deposition per track and the average ratio of all to primary tracks at the pad's location. These values, properly averaged over a region of phase space and mapped into pseudorapidity - azimuth space, constitute the double differential charged multiplicity distribution, and correcting the potentially substantial systematic errors in the result is of primary importance. The number of tracks per pad per event, called the occupancy in Phobos, is a byproduct of the reconstruction.
A subtlety arises in dealing with generated tracks that traverse more than one pad in a sensor array. These tracks must be split into segments, each lying within a single pad, so that we can determine the average energy deposited per unit track length. The layer numbers assigned to sensors in Phat are not usable in determining when a track leaves one plane, say MODC, and enters another, say MPA3, so we use a change in Sensor name, which is not a guaranteed property and which fails at some sensor joints, for this purpose. Sensor location will be used in a future version.
TPhPadMap is a data container class that provides dynamic storage and coordinate mapping functions for raw and reconstructed data. Its data member contains
TPhPadMap class methods map an array of pad data onto an array of display bins. They must be efficient because of the inherent slowness of calculating the transformations and Jacobians between Cartesian and phi-eta coordinates, which has to be done in exactly the same way for data and, with MC data input, for generated data. Therefore they are implemented as containers for both reconstructed and generated data, rather than as separate classes inheriting from a base class.
The TPhMapSet class manages an array of TPhPadMaps and has methods for accessing the individual TPhPadMaps and for creating graphs and histograms from single or combinations of TPhPadMaps.
The TPhDept class manages the data for global fits to the energy E1 deposited by a single track normally incident on a pad as a function of pad and vertex positions and total multiplicity. Its methods also fit results to obtain new values.
The TPhFSec class manages the data for global fits to the ratio of total to primary tracks as a function of pad and vertex positions and total multiplicity. Its methods also fit results to obtain new values.
The TPhMulDisp class adds axis labels, legends, and date and time stamps to plots.
The TPhOpts class provides access to ASCII files containing embedded parameter values specifying analysis options such as single track energy deposition and secondary track fractions. TPhOpts is used by TPhDept and TPhFsec.
Back to top
The most challenging part of the multiplicity reconstruction is to use front-back and azimuthal symmetries in order to reduce our dependence on models to provide the energy deposition and secondary fractions. This principle is simple: if it moves as the primary vertex changes, or as a function of azimuthal angle, it isn't physics, it's secondary background. We also use testbeam data to tune the Phobos Geant Monte Carlo parameters. However, the values required for reconstruction are averaged over the momentum and particle species spectra, and those will be known only for the averages obtained from the spectrometer. Thus the reconstruction can fail for events that are sufficiently different from the average; for example, in a plasma-enhanced model that we use for testing, a significant fraction of tracks are energetic pions at zero rapidity, and these have a very different energy deposition profile from that of a commonly-used event generator such as Hijing; using the standard profile to reconstruct some plasma events completely washes out the interesting signature of the plasma. This is an unfortunate accident, but accidents happen.
The Monte Carlo data files normally include all primary tracks but only those secondary tracks that produce hits. About 17% of charged particles in this sample produce no hits in any sensor; they include primary particles whose trajectories are outside our geometric acceptance or in inactive areas on andbetween the sensors, as well as some primaries that stop in the beam pipe or other mechanical structures. The reconstructed primary multiplicity only needs to be corrected for the small fraction that stop in the beam pipe or are past the high-eta range of our instrumentation. As we never see these tracks, the correction must be done as with any acceptance, by using simulations, with all the accompanying uncertainties. We must determine this correction and its dependence on momentum spectra and species, if significant, and assign a realistic uncertainty to it, before we can publish any results for multiplicity analysis. An additional 6% hit no multiplicity sensors, but no correction is necessary for these, since their phase space is not included in the multiplicity sensor geometry.
Back to top
h = -ln tan(q/2)
tan(q /2) = tan(q) / (1 + sqrt(1 + tan(q)2))
eh = -tan(q/2)
tan(q) = 2 tan(q/2) / (1 - tan(q/2) 2)) = 1 / sinh(h)
r = sqrt(x2 + y2)
z = r / tan(q) = r * sinh(h)
Let r be the vector to a point on a sensor plane, and rc and rv the vectors to the origin of the sensor's coordinate system and to the interaction point, respectively. Then r = r' + d, where d = rc - rv, and the local vector r' lies in the plane of the sensor. Call the unit vectors along the axes of the local coordinate system u1, u2, and u3, and choose u3 normal to the plane and pointing away from the global coordinate origin. Choose u1 to point along +z in the octagon and inward toward the beamline along the sensor's central radius in a ring. Then u2 should point along the gradient of phi in both barrel and rings, thus completing a right-handed local coordinate system for each sensor; but this simple choice has not been implemented.
A point on a sensor satisfies
(r' + d) 2 = (z' + dz) 2 / cos(q)2
where dz = zc - zv and theta is the polar angle of the point seenfrom the interaction point rv. For the octagon, the transformation equations from global eta-phi to local x'-y' coordinates are thus
y' = s * tan(f - fc)
(x' + dz) 2 / sinh(h)2 = y'2 + s2
or
x' = -dz + - s * sinh(h) / cos(f - fc)
where s == rc.u3 is the normal distance from the beamline to the sensor plane, and phic is the azimuthal angle of the center of the sensor. The hyperbola and straight line are the expected traces of a conic secton in local coordinates.
For a ring, the corresponding transform is
x'2 + y'2 = r 2
y' = r * tan(f - fc)
where r = dz * tan(q) = dz / sinh(h) and rc is taken to be on the beamline.
The basic map operation finds the fraction of (x'1..x'2, y'1..y'2) lying within (h1.. h2, f1.. f2) and converts it to the equivalent dh * df using the Jacobian del(h,f)/del(x',y'). For this we need the indefinite integral
int(dx R(x,a)) = [x*R(x,a) + - a2 * ln(x + R(x,a))]/2
where R(x,a) == sqrt(x2 + - a2)
For the octagon, the Jacobian is cos(f - fc) 3 / (s2 * cosh(h))
Back to top
The ROOT interpreter, CInt, regularly issues only two error messages. It issues the first,
*** Break *** segmentation violation
when your code or the code in your script has an error, and the second,
*** Break *** keyboard interrupt
when you press ctrl-C to try to continue. You can locate the offending interpreted code by issuing the command
.trace
before executing the script. This prints every interpreted line. If you put in printf statements as debug aids, CInt will not print anything until it interprets a newline.
The advantage of scripts is that they can be edited and re-used without leaving the root environment, either by leaving an X-window open for editing the script (this does not work well on a Mac) or by issuing the command
.!emacs -nw <script filename>
to edit the file in your command window. A variable once defined cannot be redefined, so to repair one of the most common mistakes you will have to leave ROOT, fix the error, and restart CINT.
The c++ compiler will give you compile-time errror messages. However, ROOT provides no traceback information for compiled code, and so the only way to locate the source of a crash is to use the debugger, ddd, which does not work on Abacus because of a bad installation, or to add print statements to the code. Tracing information of this sort is available from TPhMul if the relevant class global ReportLevel variables are correctly set.
If you want to know the properties of a loaded class, for example TPhMul, you can issue the command
.class TPhMul
This will print the header definition of every method in the class, with no description. You can terminate the dump by pressing ctrl-C.
The approved way of learning about Phat classes is to look in the Phat web pages, at the automatically generated document, which displays the code authors' initial comments; code developers simply look at the Phat code. Similarly, you can look here for ROOT information.
Back to top
This is done automatically by the standard Phat make; see the Phobos web page for more information.
Back to top
As of 18-Mar-1999, there are no known errors in the algorithms, and the default set of parameters, which is independent of f, seems to work well with random vertices and centralities. Work is proceeding on better tunes and the variation of parameters with different event models. More current information on status is available. Some methods described herein are not yet fully implemented. The phatloop macros used in MDC-2 run and write output datasets containing selected TPhMapSets, TGraphErrors, and Th2F's.
Back to top
One of the difficulties with CINT is that a simple enumeration, such as
enum { OK = 0, badArg = -1, unImplemented = -2 };
fails with the error messages
Error: No memory for static (null)\0\ffffffff FILE:/tmp/20970aaa LINE:1
Error: No memory for static (null)\0\ffffffff FILE:/tmp/20970aaa LINE:1
Error: No memory for static (null)\0\ffffffff FILE:/tmp/20970aaa LINE:1
*** Interpreter error recovered ***
which means that such variables have to be referred to by magic numbers, static ints, or quoted strings. There is a way to access enums that have already been defined in compiled code, by including
#pragma link C++ global <name of enumerated variable>
in create_cint for each enumerated variable; see TPhAbsRawHitFEC.h and $PHATHOME/phatdat/inc/create_cint for examples, or contact Patrick Decowski for more information.
Scheduled work:
Back to top
Appendix A: Defining the phat environment
ROOT and Phat both require several environment variables to be defined. Assuming that you are in the tcsh shell, you define these variables by using the setenv command. The fundamental variables are ROOTSYS and PHATHOME, which for the author's development version on ppcln, have the definitions
setenv ROOTSYS /home/Root/root
setenv PHATHOME ~verdier/Phat/CurrentVersion
On other systems, these directories will have different locations.
Using these definitions, you must define
setenv LD_LIBRARY_PATH $ROOTSYS/lib
setenv PHATLIB $PHATHOME/lib/`uname`
setenv PATH ${PATH}:${ROOTSYS}/bin:${PHATHOME}/bin/`uname`
where the symbol uname will be translated by the operating system's name for itself; for Linux, this is "Linux".
An alternative to redefining the PATH variable in order to be able to execute the phat module is to define phat directly using an alias:
alias phat $PHATHOME/bin/`uname`/phat
although this use of alias is strongly deprecated by Patrick Decowski.
If the setenv command is not accepted by your unix system, try typing
tcsh
to change to the tcsh shell. You can make this change permanent (recommended) by typing the change-shell command, at least on Linux:
chsh -s /bin/tcsh <username>
and following the prompts.
These variables can be defined by a setup script, if you use the slightly different form
setenv LD_LIBRARY_PATH ${ROOTSYS}/lib
setenv PHATLIB ${PHATHOME}/lib/`uname`
for these commands. If the script file is named, say, setup, you type
source setup
to make the definitions effective.
Back to top
Last edited: Tuesday, March 23, 1999 03:17:17 PM by Robin Verdier