Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:42

0001 #ifndef __HiEvtPlaneFlatten__
0002 #define __HiEvtPlaneFlatten__
0003 #include <memory>
0004 #include <iostream>
0005 #include <string>
0006 
0007 #include "FWCore/Framework/interface/Frameworkfwd.h"
0008 #include "FWCore/Framework/interface/EDProducer.h"
0009 
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/MakerMacros.h"
0012 
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 
0015 #include "DataFormats/Common/interface/Handle.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 
0018 #include "DataFormats/HeavyIonEvent/interface/EvtPlane.h"
0019 
0020 #include <vector>
0021 #include <cmath>
0022 
0023 //
0024 // class declaration
0025 //
0026 
0027 class HiEvtPlaneFlatten {
0028 public:
0029   explicit HiEvtPlaneFlatten() {
0030     hbins_ = 1;
0031     hOrder_ = 9;
0032     vorder_ = 2;  //sets default order of event plane
0033   }
0034 
0035   void init(int order,
0036             int nbins,
0037             int nvtxbins = 10,
0038             double minvtx = -25,
0039             double delvtx = 5,
0040             std::string tag = "",
0041             int vord = 2) {
0042     hOrder_ = order;  //order of flattening
0043     vorder_ = vord;   //1(v1), 2(v2), 3(v3), 4(v4)
0044     nvtxbins_ = nvtxbins;
0045     nbins_ = nbins;
0046     minvtx_ = minvtx;
0047     delvtx_ = delvtx;
0048 
0049     caloCentRefMinBin_ = -1;
0050     caloCentRefMaxBin_ = -1;
0051     hbins_ = nbins * nvtxbins_ * hOrder_;
0052     obins_ = nbins * nvtxbins_;
0053     if (hbins_ > MAXCUT) {
0054       hbins_ = 1;
0055       hOrder_ = 9;
0056     }
0057     for (int i = 0; i < hbins_; i++) {
0058       flatX_[i] = 0;
0059       flatY_[i] = 0;
0060       flatXDB_[i] = 0;
0061       flatYDB_[i] = 0;
0062       flatCnt_[i] = 0;
0063     }
0064     for (int i = 0; i < obins_; i++) {
0065       xoff_[i] = 0;
0066       yoff_[i] = 0;
0067       xoffDB_[i] = 0;
0068       yoffDB_[i] = 0;
0069       xyoffcnt_[i] = 0;
0070       xyoffmult_[i] = 0;
0071       pt_[i] = 0;
0072       pt2_[i] = 0;
0073       ptDB_[i] = 0;
0074       pt2DB_[i] = 0;
0075       ptcnt_[i] = 0;
0076     }
0077   }
0078 
0079   int cutIndx(int centbin, double vtx, int iord) const {
0080     int cut;
0081     if (centbin < 0)
0082       return -1;
0083     int ibin = centbin;
0084     int ivtx = (vtx - minvtx_) / delvtx_;
0085     if (vtx < minvtx_ || ivtx >= nvtxbins_)
0086       return -1;
0087     cut = hOrder_ * nvtxbins_ * ibin + hOrder_ * ivtx + iord;
0088     if (cut < 0 || cut >= hbins_)
0089       return -1;
0090     return cut;
0091   }
0092 
0093   int offsetIndx(int centbin, double vtx) const {
0094     int cut;
0095     if (centbin < 0)
0096       return -1;
0097     int ibin = centbin;
0098     int ivtx = (vtx - minvtx_) / delvtx_;
0099     if (ivtx < 0 || ivtx > nvtxbins_)
0100       return -1;
0101     cut = nvtxbins_ * ibin + ivtx;
0102     if (cut < 0 || cut > hbins_)
0103       return -1;
0104     return cut;
0105   }
0106 
0107   void fill(double psi, double vtx, int centbin) {
0108     if (fabs(psi) > 4)
0109       return;
0110     for (int k = 0; k < hOrder_; k++) {
0111       double fsin = sin(vorder_ * (k + 1) * psi);
0112       double fcos = cos(vorder_ * (k + 1) * psi);
0113       int indx = cutIndx(centbin, vtx, k);
0114       if (indx >= 0) {
0115         flatX_[indx] += fcos;
0116         flatY_[indx] += fsin;
0117         ++flatCnt_[indx];
0118       }
0119     }
0120   }
0121   void fillOffset(double s, double c, uint m, double vtx, int centbin) {
0122     int indx = offsetIndx(centbin, vtx);
0123     if (indx >= 0) {
0124       xoff_[indx] += c;
0125       yoff_[indx] += s;
0126       xyoffmult_[indx] += m;
0127       ++xyoffcnt_[indx];
0128     }
0129   }
0130   void fillPt(double ptval, double vtx, int centbin) {
0131     int indx = offsetIndx(centbin, vtx);
0132     if (indx >= 0) {
0133       pt_[indx] += ptval;
0134       pt2_[indx] += ptval * ptval;
0135       ++ptcnt_[indx];
0136     }
0137   }
0138 
0139   void setCaloCentRefBins(const int caloCentRefMinBin, const int caloCentRefMaxBin) {
0140     caloCentRefMinBin_ = caloCentRefMinBin;
0141     caloCentRefMaxBin_ = caloCentRefMaxBin;
0142   }
0143 
0144   double etScale(double vtx, int centbin) const {
0145     int refmin = offsetIndx(caloCentRefMinBin_, vtx);
0146     int refmax = offsetIndx(caloCentRefMaxBin_, vtx);
0147     double caloCentRefVal_ = 0;
0148     for (int i = refmin; i <= refmax; i++) {
0149       caloCentRefVal_ += getPtDB(i);
0150     }
0151     caloCentRefVal_ /= refmax - refmin + 1.;
0152     if (caloCentRefMinBin_ < 0)
0153       return 1.;
0154     int indx = offsetIndx(centbin, vtx);
0155     if (indx < 0 || caloCentRefVal_ == 0 || getPtDB(indx) == 0)
0156       return 1.;
0157     return caloCentRefVal_ / getPtDB(indx);
0158   }
0159 
0160   double getW(double pt, double vtx, int centbin) const {
0161     int indx = offsetIndx(centbin, vtx);
0162     if (indx >= 0) {
0163       double scale = etScale(vtx, centbin);
0164       double ptval = getPtDB(indx) * scale;
0165       double pt2val = getPt2DB(indx) * pow(scale, 2);
0166       double rv = pt * scale - pt2val / ptval;
0167       if (ptval > 0)
0168         return rv;
0169     }
0170     return 0.;
0171   }
0172 
0173   double getFlatPsi(double psi, double vtx, int centbin) const {
0174     double correction = 0;
0175     for (int k = 0; k < hOrder_; k++) {
0176       int indx = cutIndx(centbin, vtx, k);
0177       if (indx >= 0)
0178         correction += (2. / (double)((k + 1) * vorder_)) *
0179                       (flatXDB_[indx] * sin(vorder_ * (k + 1) * psi) - flatYDB_[indx] * cos(vorder_ * (k + 1) * psi));
0180     }
0181     psi += correction;
0182     psi = bounds(psi);
0183     psi = bounds2(psi);
0184     return psi;
0185   }
0186 
0187   double soffset(double s, double vtx, int centbin) const {
0188     int indx = offsetIndx(centbin, vtx);
0189     if (indx >= 0) {
0190       return s - yoffDB_[indx];
0191     } else {
0192       return s;
0193     }
0194   }
0195 
0196   double coffset(double c, double vtx, int centbin) const {
0197     int indx = offsetIndx(centbin, vtx);
0198     if (indx >= 0)
0199       return c - xoffDB_[indx];
0200     else
0201       return c;
0202   }
0203 
0204   double offsetPsi(double s, double c) const {
0205     double psi = atan2(s, c) / vorder_;
0206     if ((fabs(s) < 1e-4) && (fabs(c) < 1e-4))
0207       psi = 0.;
0208     psi = bounds(psi);
0209     psi = bounds2(psi);
0210     return psi;
0211   }
0212 
0213   float minCent(int indx) const {
0214     int ibin = indx / nvtxbins_;
0215     return ibin * 100. / nbins_;
0216   }
0217 
0218   float maxCent(int indx) const {
0219     int ibin = indx / nvtxbins_;
0220     return (ibin + 1) * 100. / nbins_;
0221   }
0222 
0223   double minVtx(int indx) const {
0224     int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_);
0225     return minvtx_ + ivtx * delvtx_;
0226   }
0227 
0228   double maxVtx(int indx) const {
0229     int ivtx = indx - nvtxbins_ * (int)(indx / nvtxbins_);
0230     return minvtx_ + (ivtx + 1) * delvtx_;
0231   }
0232 
0233   std::string rangeString(int indx) {
0234     char buf[120];
0235     sprintf(buf, "%5.1f < cent < %5.1f; %4.1f < vtx < %4.1f", minCent(indx), maxCent(indx), minVtx(indx), maxVtx(indx));
0236     return std::string(buf);
0237   }
0238 
0239   ~HiEvtPlaneFlatten() {}
0240   int getHBins() const { return hbins_; }
0241   int getOBins() const { return obins_; }
0242   int getNvtx() const { return nvtxbins_; }
0243   double getVtxMin() const { return minvtx_; }
0244   double getVtxMax() const { return minvtx_ + nvtxbins_ * delvtx_; }
0245   int getNcent() const { return hbins_; }
0246 
0247   double getX(unsigned int bin) const { return flatX_[bin]; }
0248   double getY(unsigned int bin) const { return flatY_[bin]; }
0249   double getXoff(unsigned int bin) const { return xoff_[bin]; }
0250   double getYoff(unsigned int bin) const { return yoff_[bin]; }
0251   double getXoffDB(unsigned int bin) const { return xoffDB_[bin]; }
0252   double getYoffDB(unsigned int bin) const { return yoffDB_[bin]; }
0253   double getXYoffcnt(unsigned int bin) const { return xyoffcnt_[bin]; }
0254   double getXYoffmult(unsigned int bin) const { return xyoffmult_[bin]; }
0255   double getPt(unsigned int bin) const { return pt_[bin]; }
0256   double getPt2(unsigned int bin) const { return pt2_[bin]; }
0257   double getPtDB(unsigned int bin) const {
0258     if (bin < MAXCUTOFF) {
0259       return ptDB_[bin];
0260     } else {
0261       return 0.;
0262     }
0263   }
0264   double getPt2DB(unsigned int bin) const {
0265     if (bin < MAXCUTOFF) {
0266       return pt2DB_[bin];
0267     } else {
0268       return 0.;
0269     }
0270   }
0271   double getPtcnt(unsigned int bin) const { return ptcnt_[bin]; }
0272   double getXDB(unsigned int bin) const { return flatXDB_[bin]; }
0273   double getYDB(unsigned int bin) const { return flatYDB_[bin]; }
0274 
0275   double getCnt(unsigned int bin) const { return flatCnt_[bin]; }
0276   void setXDB(unsigned int indx, double val) { flatXDB_[indx] = val; }
0277   void setYDB(unsigned int indx, double val) { flatYDB_[indx] = val; }
0278   void setXoffDB(unsigned int indx, double val) { xoffDB_[indx] = val; }
0279   void setYoffDB(unsigned int indx, double val) { yoffDB_[indx] = val; }
0280   void setPtDB(unsigned int indx, double val) { ptDB_[indx] = val; }
0281   void setPt2DB(unsigned int indx, double val) { pt2DB_[indx] = val; }
0282   double bounds(double ang) const {
0283     if (ang < -M_PI)
0284       ang += 2. * M_PI;
0285     if (ang > M_PI)
0286       ang -= 2. * M_PI;
0287     return ang;
0288   }
0289   double bounds2(double ang) const {
0290     double range = M_PI / (double)vorder_;
0291     while (ang < -range) {
0292       ang += 2 * range;
0293     }
0294     while (ang > range) {
0295       ang -= 2 * range;
0296     }
0297     return ang;
0298   }
0299   void setCentRes1(unsigned int bin, double res, double err) {
0300     if (bin < 100) {
0301       centRes1_[bin] = res;
0302       centResErr1_[bin] = err;
0303     }
0304   }
0305   void setCentRes2(unsigned int bin, double res, double err) {
0306     if (bin < 50) {
0307       centRes2_[bin] = res;
0308       centResErr2_[bin] = err;
0309     }
0310   }
0311   void setCentRes5(unsigned int bin, double res, double err) {
0312     if (bin < 20) {
0313       centRes5_[bin] = res;
0314       centResErr5_[bin] = err;
0315     }
0316   }
0317   void setCentRes10(unsigned int bin, double res, double err) {
0318     if (bin < 10) {
0319       centRes10_[bin] = res;
0320       centResErr10_[bin] = err;
0321     }
0322   }
0323   void setCentRes20(unsigned int bin, double res, double err) {
0324     if (bin < 5) {
0325       centRes20_[bin] = res;
0326       centResErr20_[bin] = err;
0327     }
0328   }
0329   void setCentRes25(unsigned int bin, double res, double err) {
0330     if (bin < 4) {
0331       centRes25_[bin] = res;
0332       centResErr25_[bin] = err;
0333     }
0334   }
0335   void setCentRes30(unsigned int bin, double res, double err) {
0336     if (bin < 3) {
0337       centRes30_[bin] = res;
0338       centResErr30_[bin] = err;
0339     }
0340   }
0341   void setCentRes40(unsigned int bin, double res, double err) {
0342     if (bin < 2) {
0343       centRes40_[bin] = res;
0344       centResErr40_[bin] = err;
0345     }
0346   }
0347 
0348   double getCentRes1(unsigned int bin) const {
0349     if (bin < 100) {
0350       return centRes1_[bin];
0351     } else {
0352       return 0.;
0353     }
0354   }
0355   double getCentRes2(unsigned int bin) const {
0356     if (bin < 50) {
0357       return centRes2_[bin];
0358     } else {
0359       return 0.;
0360     }
0361   }
0362   double getCentRes5(unsigned int bin) const {
0363     if (bin < 20) {
0364       return centRes5_[bin];
0365     } else {
0366       return 0.;
0367     }
0368   }
0369   double getCentRes10(unsigned int bin) const {
0370     if (bin < 10) {
0371       return centRes10_[bin];
0372     } else {
0373       return 0.;
0374     }
0375   }
0376   double getCentRes20(unsigned int bin) const {
0377     if (bin < 5) {
0378       return centRes20_[bin];
0379     } else {
0380       return 0.;
0381     }
0382   }
0383   double getCentRes25(unsigned int bin) const {
0384     if (bin < 4) {
0385       return centRes25_[bin];
0386     } else {
0387       return 0.;
0388     }
0389   }
0390   double getCentRes30(unsigned int bin) const {
0391     if (bin < 3) {
0392       return centRes30_[bin];
0393     } else {
0394       return 0.;
0395     }
0396   }
0397   double getCentRes40(unsigned int bin) const {
0398     if (bin < 2) {
0399       return centRes40_[bin];
0400     } else {
0401       return 0.;
0402     }
0403   }
0404 
0405   double getCentResErr1(unsigned int bin) const {
0406     if (bin < 100) {
0407       return centResErr1_[bin];
0408     } else {
0409       return 0.;
0410     }
0411   }
0412   double getCentResErr2(unsigned int bin) const {
0413     if (bin < 50) {
0414       return centResErr2_[bin];
0415     } else {
0416       return 0.;
0417     }
0418   }
0419   double getCentResErr5(unsigned int bin) const {
0420     if (bin < 20) {
0421       return centResErr5_[bin];
0422     } else {
0423       return 0.;
0424     }
0425   }
0426   double getCentResErr10(unsigned int bin) const {
0427     if (bin < 10) {
0428       return centResErr10_[bin];
0429     } else {
0430       return 0.;
0431     }
0432   }
0433   double getCentResErr20(unsigned int bin) const {
0434     if (bin < 5) {
0435       return centResErr20_[bin];
0436     } else {
0437       return 0.;
0438     }
0439   }
0440   double getCentResErr25(unsigned int bin) const {
0441     if (bin < 4) {
0442       return centResErr25_[bin];
0443     } else {
0444       return 0.;
0445     }
0446   }
0447   double getCentResErr30(unsigned int bin) const {
0448     if (bin < 3) {
0449       return centResErr30_[bin];
0450     } else {
0451       return 0.;
0452     }
0453   }
0454   double getCentResErr40(unsigned int bin) const {
0455     if (bin < 2) {
0456       return centResErr40_[bin];
0457     } else {
0458       return 0.;
0459     }
0460   }
0461 
0462 private:
0463   static const int MAXCUT = 10000;
0464   static const int MAXCUTOFF = 1000;
0465 
0466   double flatX_[MAXCUT];
0467   double flatY_[MAXCUT];
0468   double flatXDB_[MAXCUT];
0469   double flatYDB_[MAXCUT];
0470   double flatCnt_[MAXCUT];
0471 
0472   double xoff_[MAXCUTOFF];
0473   double yoff_[MAXCUTOFF];
0474   double xoffDB_[MAXCUTOFF];
0475   double yoffDB_[MAXCUTOFF];
0476   double xyoffcnt_[MAXCUTOFF];
0477   uint xyoffmult_[MAXCUTOFF];
0478 
0479   double pt_[MAXCUTOFF];
0480   double pt2_[MAXCUTOFF];
0481   double ptDB_[MAXCUTOFF];
0482   double pt2DB_[MAXCUTOFF];
0483   double ptcnt_[MAXCUTOFF];
0484 
0485   double centRes1_[100];
0486   double centResErr1_[100];
0487 
0488   double centRes2_[50];
0489   double centResErr2_[50];
0490 
0491   double centRes5_[20];
0492   double centResErr5_[20];
0493 
0494   double centRes10_[10];
0495   double centResErr10_[10];
0496 
0497   double centRes20_[5];
0498   double centResErr20_[5];
0499 
0500   double centRes25_[4];
0501   double centResErr25_[4];
0502 
0503   double centRes30_[3];
0504   double centResErr30_[3];
0505 
0506   double centRes40_[2];
0507   double centResErr40_[2];
0508 
0509   int nvtxbins_;
0510   int nbins_;
0511   double minvtx_;
0512   double delvtx_;
0513   int hOrder_;             //flattening order
0514   int hbins_;              //number of bins needed for flattening
0515   int obins_;              //number of (x,y) offset bins
0516   int vorder_;             //order of flattened event plane
0517   int caloCentRefMinBin_;  //min ref centrality bin for calo weight scale
0518   int caloCentRefMaxBin_;  //max ref centrality bin for calo weight scale
0519 };
0520 
0521 #endif