Home » Blogs » lisa's blog

Using the EPD event-plane finder

Updated on Wed, 2018-10-24 15:56. Originally created by lisa on 2018-10-16 10:21.

How to find the Event Plane with the STAR EPD code StEpdEpFinder

Contents:

General points

  • All code described here is in STAR cvs under StRoot/StEpdUtil/
  • The object that you instantiate and use in your code is StEpdEpFinder, which returns the results in a StEpdEpInfo object.  This object contains all angles, Q-vectors, ring-by-ring information, everything.
  • The maximum "order" of the event plane calculated is set by a pre-compiler directive in the header file StEpdEpInfo.h.  The minimum order is always 1.
    • #define _EpOrderMax 3    // EP is calculated to order 3
  • Whenever the user specifies the "order" of the event plane
    • order=1 means first-order plane (i.e. "directed flow" plane)
    • order=2 means second-order plane (i.e. "elliptic flow" plane)
    • etc
    • There is no "order=0" plane.  Internally in the code, the arrays start at zero as you might expect, but all of this is hidden from the user.
  • The StEpdEpFinder applies corrections (inefficiency weighting, angle shifting, etc) to the raw data; it reads these in from a file.  The exact same code also generates these corrections and stores them into a file.  You don't really have to do anything.  The code will generate the corrections and put them into a file whether you ask it to, or not.  All you have to do, is make sure you give the code the proper input.  This iterative procedure is common in many event-plane finder codes.  Details are discussed below.
  • Contact Mike Lisa (lisa.1@osu.edu) with problems.

Before you run:

  • Decide the maximum order that you want to calculate (n=3 is the "triangular flow" plane) and set the compiler variable _EpOrderMax in StEpdEpInfo.h.  Compile the code.
  • Decide how many centrality bins you will use, and your scheme for identifying them using an integer.  (If you don't care about centrality, you can just have one big bin, and every event has e.g. centId=1.)

Settings that you make while and just after instantiating the StEpdEpFinder in your code

  • Specify the number of centrality bins and filename for the correction information
    StEpdEpFinder::StEpdEpFinder(int nCentralityBins=10, const char* CorrectionFileName="StEpdEpFinderCorrectionHistograms_INPUT.root");
    • If you don't yet have a file with correction factors, don't worry.  The code will simply not use correction factors, but will work fine.
  • Select your data source.  The code will be getting data from the StEpdHit (or StMuEpdHit or StPicoEpdHit) objects.
    • StEpdEpFinder::SetEpdHitFormat(int format);    // format=0/1/2 for StEpdHit/StMuEpdHit/StPicoEpdHit
  • Select weighting parameters (optional but recommended).   As discussed extensively in the EPD group (e.g. here), it is best to weight an EPD hit by applying a threshold and a maximum (ceiling) value to the nMIP data.
    • StEpdEpFinder::SetnMipThreshold(double thresh);    // recommended value 0.3
    • StEpdEpFinder::SetMaxTileWeight(double MAX);        // recommended ~1 for low multiplicity (e.g. BES) to 3 for high multiplicity (e.g. 200 GeV)
    • you may wish to play around with these settings, but we have found the dependence of the EP resolution to be mild
  • Set ring or eta weights (optional but can be very important).   Some regions of phasespace are more sensitive to the event-plane than others, so should be "weighted" more highly.  Indeed, some regions should get a negative weight!  To handle this, you can give the code weights based either on ring (1..16) or eta range.  Eta range is recommended, and is probably only really important for first-order EP.  The weights will depend on EP order (and on centrality, if you use eta weighting).
    • To use weighting by ring:
      StEpdEpFinder::SetRingWeights(int order, double* RingWeights);    // RingWeights is a 1D array of 16 elements.
    • To use weighting according to eta:
      StEpdEpFinder::SetEtaWeights(int order, TH2D EtaWeight);   // histogram is binned in |eta| and centrality, according to the the centrality scheme you are using.  This way, you choose the centrality binning and the eta binning you want.

Getting the planes

For every event, just call StEpdEpFinder::Results(TClonesArray* EpdHits, TVector3 primVertex, int CentralityID);
  • EpdHits contains the StEpdHit/StMuEpdHit/StPicoEpdHit objects, as specified by SetEpdHitFormat.  (See above).
  • The primary vertex is required, to determine the pseudorapidity (and potentially azimuthal angle) of particles hitting the tile.  A straight-line approximation is used.
  • The CentralityID, in whatever scheme you have chosen (remember?), is used both for applying the proper correction factors, and for incrementing the proper correction factors.
This method returns a StEpdEpInfo object, which contains all the information you could ever want.  (Reason: it is easy to just calculate all this stuff at the same time.  Just giving the user the entire bundle all at once, avoids unnecessary re-calculation if he asks for more than one item of information.)  There is EP information for the full EPD, for East/West only, and for rings.  You have Q-vectors and EP angles.  You have this information for various levels of correction.  Probably best if you explore the information yourself: see the Doxygen page. 
  • If you just want "the best EP angle, fully corrected and using all of the information," then you want
    double StEpdEpInfo::FullPhiWeightedAndShiftedPsi(int order)

Corrections

The code does two corrections.
  1. The "Phi-Weight" correction, which simply suppresses "hot" and enhances "cold" tiles, forcing all tiles in a given ring to have the same average signal.
  2. The "Psi-shifting" correction, which implements equation (6) of this famous paper by Poskanzer and Voloshin
The factors for these two corrections are contained in histograms in a TFile that StEpdEpFinder reads in upon initial configuration.  You give the TFile name when instantiating StEpdEpFinder (in the ctor).  The correction factors themselves are in fact calculated by the StEpdEpFinder, automatically.  So, the way you run is:
  1. Run over the data once.  In addition to whatever output you produce, the code will create a file StEpdEpFinderCorrectionHistograms_OUTPUT.root
  2. mv StEpdEpFinderCorrectionHistograms_OUTPUT.root StEpdEpFinderCorrectionHistograms_INPUT.root
    • At this point, the Phi-Weight corrections in the file are correct, but not the Psi-shift corrections
  3. Run over the data again, exactly as the first time.
  4. And again mv StEpdEpFinderCorrectionHistograms_OUTPUT.root StEpdEpFinderCorrectionHistograms_INPUT.root
    • At this point, both the Phi-Weight and the Psi-shift corrections are correct.
  5. And now you are good to go.  No need to rename the file anymore.  Just run as many times as you want.
As it turns out, these corrections are usually quite mild in the EPD, which comes pre-calibrated more or less "out of the box," is pretty darned symmetric, and does not have hot and cold spots.  Unlike (ahem) some detectors...

Nuts and bolts 1 - code snips (obviously, just a specific example)

This is just some relevant snips from a very stripped down analysis code, which you may download-- without warranty and in as-is, undocumented condition-- here.

Put this in your header...
#include "StEpdUtil/StEpdEpFinder.h"

Before doing event-by-event processing...


1) When you have the picoDst, get a pointer to the TClonesArray of StPicoEpdHit objects.  (Or StMuEpdHit objects, if usig muDsts.)  You will be passing this to the EP finder.
  mEpdHits = new TClonesArray("StPicoEpdHit");
  mPicoDst->SetBranchStatus("EpdHit*",1,&found);   // note you need the asterisk                                                                                
  cout << "EpdHit Branch returned found= " << found << endl;
  mPicoDst->SetBranchAddress("EpdHit",&mEpdHits);


2) Instantiate and configure the EP finder...
  StEpdEpFinder* mEpFinder = new StEpdEpFinder();
  mEpFinder->SetnMipThreshold(0.3);    // recommended by EPD group
  mEpFinder->SetMaxTileWeight(2.0);     // recommended by EPD group
  mEpFinder->SetEpdHitFormat(2);         // 2=pico


3) For first-order flow using the 27 GeV Au+Au data, I use eta-based weights as follows: (c.f. here)
  double lin[9] = {-1.950, -1.900, -1.850, -1.706, -1.438, -1.340, -1.045, -0.717, -0.700};
  double cub[9] = {0.1608, 0.1600, 0.1600, 0.1595, 0.1457, 0.1369, 0.1092, 0.0772, 0.0700};
  TH2D wt("Order1etaWeight","Order1etaWeight",100,1.5,6.5,10,0,10);
  for (int ix=1; ix<101; ix++){
    for (int iy=1; iy<11; iy++){
      double eta = wt.GetXaxis()->GetBinCenter(ix);
      wt.SetBinContent(ix,iy,lin[iy-1]*eta+cub[iy-1]*pow(eta,3));
    }
  }
  mEpFinder->SetEtaWeights(1,wt);

Then, in your event loop,

just...
  StEpdEpInfo result = mEpFinder->Results(mEpdHits,PV,CentId);

...and do with that information what you will.  E.g. fill some histograms like this:
  for (int order=1; order<4; order++){
    mEpdAveCos[order-1][0]->Fill(double(CentId),cos((double)order*(result.EastRawPsi(order)-result.WestRawPsi(order))));
    mEpdAveCos[order-1][1]->Fill(double(CentId),cos((double)order*(result.EastPhiWeightedAndShiftedPsi(order)-result.WestPhiWeightedAndShiftedPsi(order))));
    mEpdEwPsi[order-1][0]->Fill(result.EastRawPsi(order),result.WestRawPsi(order));
    mEpdEwPsi[order-1][1]->Fill(result.EastPhiWeightedAndShiftedPsi(order),result.WestPhiWeightedAndShiftedPsi(order));
  }

Nuts and bolts 2 - iterative running

At the command line, I do this:
first iteration:
./AnalyzePico 2000 /Volumes/MikeLisaPassport/NewPicos/Au27Au2018picos/27GeV_production/AllDays.list 0
mv StEpdEpFinderCorrectionHistograms_OUTPUT.root StEpdEpFinderCorrectionHistograms_INPUT.root
second iteration:
./AnalyzePico 2000 /Volumes/MikeLisaPassport/NewPicos/Au27Au2018picos/27GeV_production/AllDays.list 0
mv StEpdEpFinderCorrectionHistograms_OUTPUT.root StEpdEpFinderCorrectionHistograms_INPUT.root

Now, done.  Run some more, change your code to do other things with the event plane, whatever.  No need to overwrite your INPUT correction file anymore, unless
1) you change the dataset you are analyzing
2) you want to change the parameters (e.g. EPD threshold or weighting scheme)

./AnalyzePico 2000 /Volumes/MikeLisaPassport/NewPicos/Au27Au2018picos/27GeV_production/AllDays.list 0



»