Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:18 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!