HiEvtPlaneFlatten.h

Go to the documentation of this file.

#ifndef __HiEvtPlaneFlatten__
#define __HiEvtPlaneFlatten__
#include <memory>
#include <iostream>
#include <string>

#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDProducer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"

#include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"

#include <vector>
#include <cmath>

//
// class declaration
//

class HiEvtPlaneFlatten {
public:
  explicit HiEvtPlaneFlatten() {
    hbins_ = 1;
    hOrder_ = 9;
    vorder_ = 2;  //sets default order of event plane
  }

  void init(int order,
            int nbins,
            int nvtxbins = 10,
            double minvtx = -25,
            double delvtx = 5,
            std::string tag = "",
            int vord = 2) {
    hOrder_ = order;  //order of flattening
    vorder_ = vord;   //1(v1), 2(v2), 3(v3), 4(v4)
    nvtxbins_ = nvtxbins;
    nbins_ = nbins;
    minvtx_ = minvtx;
    delvtx_ = delvtx;

    caloCentRefMinBin_ = -1;
    caloCentRefMaxBin_ = -1;
    hbins_ = nbins * nvtxbins_ * hOrder_;
    obins_ = nbins * nvtxbins_;
    if (hbins_ > MAXCUT) {
      hbins_ = 1;
      hOrder_ = 9;
    }
    for (int i = 0; i < hbins_; i++) {
      flatX_[i] = 0;
      flatY_[i] = 0;
      flatXDB_[i] = 0;
      flatYDB_[i] = 0;
      flatCnt_[i] = 0;
    }
    for (int i = 0; i < obins_; i++) {
      xoff_[i] = 0;
      yoff_[i] = 0;
      xoffDB_[i] = 0;
      yoffDB_[i] = 0;
      xyoffcnt_[i] = 0;
      xyoffmult_[i] = 0;
      pt_[i] = 0;
      pt2_[i] = 0;
      ptDB_[i] = 0;
      pt2DB_[i] = 0;
      ptcnt_[i] = 0;
    }
  }

  int cutIndx(int centbin, double vtx, int iord) const {
    int cut;
    if (centbin < 0)
      return -1;
    int ibin = centbin;
    int ivtx = (vtx - minvtx_) / delvtx_;
    if (vtx < minvtx_ || ivtx >= nvtxbins_)
      return -1;
    cut = hOrder_ * nvtxbins_ * ibin + hOrder_ * ivtx + iord;
    if (cut < 0 || cut >= hbins_)
      return -1;
    return cut;
  }

  int offsetIndx(int centbin, double vtx) const {
    int cut;
    if (centbin < 0)
      return -1;
    int ibin = centbin;
    int ivtx = (vtx - minvtx_) / delvtx_;
    if (ivtx < 0 || ivtx > nvtxbins_)
      return -1;
    cut = nvtxbins_ * ibin + ivtx;
    if (cut < 0 || cut > hbins_)
      return -1;
    return cut;
  }

  void fill(double psi, double vtx, int centbin) {
    if (fabs(psi) > 4)
      return;
    for (int k = 0; k < hOrder_; k++) {
      double fsin = sin(vorder_ * (k + 1) * psi);
      double fcos = cos(vorder_ * (k + 1) * psi);
      int indx = cutIndx(centbin, vtx, k);
      if (indx >= 0) {
        flatX_[indx] += fcos;
        flatY_[indx] += fsin;
        ++flatCnt_[indx];
      }
    }
  }
  void fillOffset(double s, double c, uint m, double vtx, int centbin) {
    int indx = offsetIndx(centbin, vtx);
    if (indx >= 0) {
      xoff_[indx] += c;
      yoff_[indx] += s;
      xyoffmult_[indx] += m;
      ++xyoffcnt_[indx];
    }
  }
  void fillPt(double ptval, double vtx, int centbin) {
    int indx = offsetIndx(centbin, vtx);
    if (indx >= 0) {
      pt_[indx] += ptval;
      pt2_[indx] += ptval * ptval;
      ++ptcnt_[indx];
    }
  }

  void setCaloCentRefBins(const int caloCentRefMinBin, const int caloCentRefMaxBin) {
    caloCentRefMinBin_ = caloCentRefMinBin;
    caloCentRefMaxBin_ = caloCentRefMaxBin;
  }

  double etScale(double vtx, int centbin) const {
    int refmin = offsetIndx(caloCentRefMinBin_, vtx);
    int refmax = offsetIndx(caloCentRefMaxBin_, vtx);
    double caloCentRefVal_ = 0;
    for (int i = refmin; i <= refmax; i++) {
      caloCentRefVal_ += getPtDB(i);
    }
    caloCentRefVal_ /= refmax - refmin + 1.;
    if (caloCentRefMinBin_ < 0)
      return 1.;
    int indx = offsetIndx(centbin, vtx);
    if (indx < 0 || caloCentRefVal_ == 0 || getPtDB(indx) == 0)
      return 1.;
    return caloCentRefVal_ / getPtDB(indx);
  }

  double getW(double pt, double vtx, int centbin) const {
    int indx = offsetIndx(centbin, vtx);
    if (indx >= 0) {
      double scale = etScale(vtx, centbin);
      double ptval = getPtDB(indx) * scale;
      double pt2val = getPt2DB(indx) * pow(scale, 2);
      double rv = pt * scale - pt2val / ptval;
      if (ptval > 0)
        return rv;
    }
    return 0.;
  }

  double getFlatPsi(double psi, double vtx, int centbin) const {
    double correction = 0;
    for (int k = 0; k < hOrder_; k++) {
      int indx = cutIndx(centbin, vtx, k);
      if (indx >= 0)
        correction += (2. / (double)((k + 1) * vorder_)) *
                      (flatXDB_[indx] * sin(vorder_ * (k + 1) * psi) - flatYDB_[indx] * cos(vorder_ * (k + 1) * psi));
    }
    psi += correction;
    psi = bounds(psi);
    psi = bounds2(psi);
    return psi;
  }

  double soffset(double s, double vtx, int centbin) const {
    int indx = offsetIndx(centbin, vtx);
    if (indx >= 0) {
      return s - yoffDB_[indx];
    } else {
      return s;
    }
  }

  double coffset(double c, double vtx, int centbin) const {
    int indx = offsetIndx(centbin, vtx);
    if (indx >= 0)
      return c - xoffDB_[indx];
    else
      return c;
  }

  double offsetPsi(double s, double c) const {
    double psi = atan2(s, c) / vorder_;
    if ((fabs(s) < 1e-4) && (fabs(c) < 1e-4))
      psi = 0.;
    psi = bounds(psi);
    psi = bounds2(psi);
    return psi;
  }

  float minCent(int indx) const {
    int ibin = indx / nvtxbins_;
    return ibin * 100. / nbins_;
  }

  float maxCent(int indx) const {
    int ibin = indx / nvtxbins_;
    return (ibin + 1) * 100. / nbins_;
  }

  double minVtx(int indx) const {
    int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_);
    return minvtx_ + ivtx * delvtx_;
  }

  double maxVtx(int indx) const {
    int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_);
    return minvtx_ + (ivtx + 1) * delvtx_;
  }

  std::string rangeString(int indx) {
    char buf[120];
    sprintf(buf, "%5.1f < cent < %5.1f; %4.1f < vtx < %4.1f", minCent(indx), maxCent(indx), minVtx(indx), maxVtx(indx));
    return std::string(buf);
  }

  ~HiEvtPlaneFlatten() {}
  int getHBins() const { return hbins_; }
  int getOBins() const { return obins_; }
  int getNvtx() const { return nvtxbins_; }
  double getVtxMin() const { return minvtx_; }
  double getVtxMax() const { return minvtx_ + nvtxbins_ * delvtx_; }
  int getNcent() const { return hbins_; }

  double getX(unsigned int bin) const { return flatX_[bin]; }
  double getY(unsigned int bin) const { return flatY_[bin]; }
  double getXoff(unsigned int bin) const { return xoff_[bin]; }
  double getYoff(unsigned int bin) const { return yoff_[bin]; }
  double getXoffDB(unsigned int bin) const { return xoffDB_[bin]; }
  double getYoffDB(unsigned int bin) const { return yoffDB_[bin]; }
  double getXYoffcnt(unsigned int bin) const { return xyoffcnt_[bin]; }
  double getXYoffmult(unsigned int bin) const { return xyoffmult_[bin]; }
  double getPt(unsigned int bin) const { return pt_[bin]; }
  double getPt2(unsigned int bin) const { return pt2_[bin]; }
  double getPtDB(unsigned int bin) const {
    if (bin < MAXCUTOFF) {
      return ptDB_[bin];
    } else {
      return 0.;
    }
  }
  double getPt2DB(unsigned int bin) const {
    if (bin < MAXCUTOFF) {
      return pt2DB_[bin];
    } else {
      return 0.;
    }
  }
  double getPtcnt(unsigned int bin) const { return ptcnt_[bin]; }
  double getXDB(unsigned int bin) const { return flatXDB_[bin]; }
  double getYDB(unsigned int bin) const { return flatYDB_[bin]; }

  double getCnt(unsigned int bin) const { return flatCnt_[bin]; }
  void setXDB(unsigned int indx, double val) { flatXDB_[indx] = val; }
  void setYDB(unsigned int indx, double val) { flatYDB_[indx] = val; }
  void setXoffDB(unsigned int indx, double val) { xoffDB_[indx] = val; }
  void setYoffDB(unsigned int indx, double val) { yoffDB_[indx] = val; }
  void setPtDB(unsigned int indx, double val) { ptDB_[indx] = val; }
  void setPt2DB(unsigned int indx, double val) { pt2DB_[indx] = val; }
  double bounds(double ang) const {
    if (ang < -M_PI)
      ang += 2. * M_PI;
    if (ang > M_PI)
      ang -= 2. * M_PI;
    return ang;
  }
  double bounds2(double ang) const {
    double range = M_PI / (double)vorder_;
    while (ang < -range) {
      ang += 2 * range;
    }
    while (ang > range) {
      ang -= 2 * range;
    }
    return ang;
  }
  void setCentRes1(unsigned int bin, double res, double err) {
    if (bin < 100) {
      centRes1_[bin] = res;
      centResErr1_[bin] = err;
    }
  }
  void setCentRes2(unsigned int bin, double res, double err) {
    if (bin < 50) {
      centRes2_[bin] = res;
      centResErr2_[bin] = err;
    }
  }
  void setCentRes5(unsigned int bin, double res, double err) {
    if (bin < 20) {
      centRes5_[bin] = res;
      centResErr5_[bin] = err;
    }
  }
  void setCentRes10(unsigned int bin, double res, double err) {
    if (bin < 10) {
      centRes10_[bin] = res;
      centResErr10_[bin] = err;
    }
  }
  void setCentRes20(unsigned int bin, double res, double err) {
    if (bin < 5) {
      centRes20_[bin] = res;
      centResErr20_[bin] = err;
    }
  }
  void setCentRes25(unsigned int bin, double res, double err) {
    if (bin < 4) {
      centRes25_[bin] = res;
      centResErr25_[bin] = err;
    }
  }
  void setCentRes30(unsigned int bin, double res, double err) {
    if (bin < 3) {
      centRes30_[bin] = res;
      centResErr30_[bin] = err;
    }
  }
  void setCentRes40(unsigned int bin, double res, double err) {
    if (bin < 2) {
      centRes40_[bin] = res;
      centResErr40_[bin] = err;
    }
  }

  double getCentRes1(unsigned int bin) const {
    if (bin < 100) {
      return centRes1_[bin];
    } else {
      return 0.;
    }
  }
  double getCentRes2(unsigned int bin) const {
    if (bin < 50) {
      return centRes2_[bin];
    } else {
      return 0.;
    }
  }
  double getCentRes5(unsigned int bin) const {
    if (bin < 20) {
      return centRes5_[bin];
    } else {
      return 0.;
    }
  }
  double getCentRes10(unsigned int bin) const {
    if (bin < 10) {
      return centRes10_[bin];
    } else {
      return 0.;
    }
  }
  double getCentRes20(unsigned int bin) const {
    if (bin < 5) {
      return centRes20_[bin];
    } else {
      return 0.;
    }
  }
  double getCentRes25(unsigned int bin) const {
    if (bin < 4) {
      return centRes25_[bin];
    } else {
      return 0.;
    }
  }
  double getCentRes30(unsigned int bin) const {
    if (bin < 3) {
      return centRes30_[bin];
    } else {
      return 0.;
    }
  }
  double getCentRes40(unsigned int bin) const {
    if (bin < 2) {
      return centRes40_[bin];
    } else {
      return 0.;
    }
  }

  double getCentResErr1(unsigned int bin) const {
    if (bin < 100) {
      return centResErr1_[bin];
    } else {
      return 0.;
    }
  }
  double getCentResErr2(unsigned int bin) const {
    if (bin < 50) {
      return centResErr2_[bin];
    } else {
      return 0.;
    }
  }
  double getCentResErr5(unsigned int bin) const {
    if (bin < 20) {
      return centResErr5_[bin];
    } else {
      return 0.;
    }
  }
  double getCentResErr10(unsigned int bin) const {
    if (bin < 10) {
      return centResErr10_[bin];
    } else {
      return 0.;
    }
  }
  double getCentResErr20(unsigned int bin) const {
    if (bin < 5) {
      return centResErr20_[bin];
    } else {
      return 0.;
    }
  }
  double getCentResErr25(unsigned int bin) const {
    if (bin < 4) {
      return centResErr25_[bin];
    } else {
      return 0.;
    }
  }
  double getCentResErr30(unsigned int bin) const {
    if (bin < 3) {
      return centResErr30_[bin];
    } else {
      return 0.;
    }
  }
  double getCentResErr40(unsigned int bin) const {
    if (bin < 2) {
      return centResErr40_[bin];
    } else {
      return 0.;
    }
  }

private:
  static const int MAXCUT = 10000;
  static const int MAXCUTOFF = 1000;

  double flatX_[MAXCUT];
  double flatY_[MAXCUT];
  double flatXDB_[MAXCUT];
  double flatYDB_[MAXCUT];
  double flatCnt_[MAXCUT];

  double xoff_[MAXCUTOFF];
  double yoff_[MAXCUTOFF];
  double xoffDB_[MAXCUTOFF];
  double yoffDB_[MAXCUTOFF];
  double xyoffcnt_[MAXCUTOFF];
  uint xyoffmult_[MAXCUTOFF];

  double pt_[MAXCUTOFF];
  double pt2_[MAXCUTOFF];
  double ptDB_[MAXCUTOFF];
  double pt2DB_[MAXCUTOFF];
  double ptcnt_[MAXCUTOFF];

  double centRes1_[100];
  double centResErr1_[100];

  double centRes2_[50];
  double centResErr2_[50];

  double centRes5_[20];
  double centResErr5_[20];

  double centRes10_[10];
  double centResErr10_[10];

  double centRes20_[5];
  double centResErr20_[5];

  double centRes25_[4];
  double centResErr25_[4];

  double centRes30_[3];
  double centResErr30_[3];

  double centRes40_[2];
  double centResErr40_[2];

  int nvtxbins_;
  int nbins_;
  double minvtx_;
  double delvtx_;
  int hOrder_;             //flattening order
  int hbins_;              //number of bins needed for flattening
  int obins_;              //number of (x,y) offset bins
  int vorder_;             //order of flattened event plane
  int caloCentRefMinBin_;  //min ref centrality bin for calo weight scale
  int caloCentRefMaxBin_;  //max ref centrality bin for calo weight scale
};

#endif