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
0016
0017
0018 class HiEvtPlaneFlatten {
0019 public:
0020 explicit HiEvtPlaneFlatten() {
0021 hbins_ = 1;
0022 hOrder_ = 9;
0023 vorder_ = 2;
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;
0034 vorder_ = vord;
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_;
0505 int hbins_;
0506 int obins_;
0507 int vorder_;
0508 int caloCentRefMinBin_;
0509 int caloCentRefMaxBin_;
0510 };
0511
0512 #endif