EcalLaserAnalyzer/src/ME.cc

0001 #include <cassert>
0002 #include <cstdio>
0003 #include <cstdlib>
0004 #include <iostream>
0005 #include <ctime>
0006
0007 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/ME.h"
0008 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEGeom.h"
0009 #include "CalibCalorimetry/EcalLaserAnalyzer/interface/MEChannel.h"
0010
0011 //GHM ClassImp(ME)
0012
0013 using namespace std;
0014
0015 const ME::TimeStamp ME::kLowMask = 0xFFFFFFFF;
0016
0017 const TString ME::granularity[ME::iSizeG] = {"R", "SM", "LMR", "LMM", "SC", "C"};
0018
0019 const TString ME::APDPrimVar[ME::iSizeAPD] = {"FLAG",
0020                                               "MEAN",
0021                                               "RMS",
0022                                               "M3",
0023                                               "APD_OVER_PNA_MEAN",
0024                                               "APD_OVER_PNA_RMS",
0025                                               "APD_OVER_PNA_M3",
0026                                               "APD_OVER_PNB_MEAN",
0027                                               "APD_OVER_PNB_RMS",
0028                                               "APD_OVER_PNB_M3",
0029                                               "APD_OVER_PN_MEAN",
0030                                               "APD_OVER_PN_RMS",
0031                                               "APD_OVER_PN_M3",
0032                                               "SHAPE_COR",
0033                                               "ALPHA",
0034                                               "BETA",
0035                                               "TIME_MEAN",
0036                                               "TIME_RMS",
0037                                               "TIME_M3",
0038                                               "TIME_NEVT"};
0039
0040 const TString ME::PNPrimVar[ME::iSizePN] = {
0041     "FLAG",
0042     "MEAN",
0043     "RMS",
0044     "M3",
0045     "PNA_OVER_PNB_MEAN",
0046     "PNA_OVER_PNB_RMS",
0047     "PNA_OVER_PNB_M3",
0048 };
0049
0050 const TString ME::MTQPrimVar[ME::iSizeMTQ] = {
0051     "FIT_METHOD", "MTQ_AMPL", "MTQ_TIME", "MTQ_RISE", "MTQ_FWHM", "MTQ_FW20", "MTQ_FW80", "MTQ_SLIDING"};
0052
0053 const TString ME::TPAPDPrimVar[ME::iSizeTPAPD] = {"FLAG", "MEAN", "RMS", "M3", "NEVT"};
0054
0055 const TString ME::TPPNPrimVar[ME::iSizeTPPN] = {"GAIN", "MEAN", "RMS", "M3"};
0056
0057 const TString ME::type[ME::iSizeT] = {"Laser", "TestPulse"};
0058
0059 const TString ME::color[ME::iSizeC] = {"Blue", "Green", "Red", "IRed", "LED1", "LED2"};
0060
0061 std::vector<MEChannel*> ME::_trees = std::vector<MEChannel*>(4, (MEChannel*)nullptr);
0062
0063 const bool ME::useElectronicNumbering = false;
0064
0065 TString ME::lmdataPath(int lmr) {
0066   TString out_(std::getenv("MELMDAT"));
0067   out_ += "/";
0068   out_ += ME::smName(lmr);
0069   out_ += "/";
0070   return out_;
0071 }
0072
0073 TString ME::primPath(int lmr) {
0074   TString out_(std::getenv("MESTORE"));
0075   out_ += "/";
0076   out_ += ME::smName(lmr);
0077   out_ += "/";
0078   return out_;
0079 }
0080
0081 TString ME::path() { return TString(std::getenv("MUSECAL")) + "/"; }
0082
0083 TString ME::rootFileName(ME::Header header, ME::Settings settings) {
0084   // get the laser monitoring region and super-module
0085   int lmr_ = ME::lmr(header.dcc, header.side);
0086   TString outfile_ = primPath(lmr_);
0087   outfile_ += "LMF_";
0088   outfile_ += ME::smName(lmr_);
0089   outfile_ += "_";
0090   outfile_ += header.side;
0091   if (settings.type == ME::iLaser) {
0092     switch (settings.wavelength) {
0093       case iBlue:
0094         outfile_ += "_BlueLaser";
0095         break;
0096       case iGreen:
0097         outfile_ += "_GreenLaser";
0098         break;
0099       case iRed:
0100         outfile_ += "_RedLaser";
0101         break;
0102       case iIRed:
0103         outfile_ += "_IRedLaser";
0104         break;
0105       default:
0106         break;
0107     }
0108   } else if (settings.type == ME::iTestPulse) {
0109     outfile_ += "_testPulse";
0110   }
0111   outfile_ += "_";
0112   outfile_ += header.rundir.c_str();
0113   outfile_ += "_TS";
0114   outfile_ += header.ts_beg;
0115   outfile_ += ".root";
0116   return outfile_;
0117 }
0118
0119 TString ME::runListName(int lmr, int type, int color) {
0120   TString outfile_ = primPath(lmr);
0121   if (type == iLaser) {
0122     outfile_ += "runlist_";
0123     switch (color) {
0124       case ME::iBlue:
0125         outfile_ += "Blue_";
0126         break;
0127       case ME::iGreen:
0128         outfile_ += "Red_";
0129         break;
0130       case ME::iRed:
0131         outfile_ += "Red_";
0132         break;
0133       case ME::iIRed:
0134         outfile_ += "IRed_";
0135         break;
0136       default:
0137         abort();
0138     }
0139     outfile_ += "Laser";
0140   } else if (type == iTestPulse) {
0141     outfile_ += "runlist_Test_Pulse";
0142   }
0143   return outfile_;
0144 }
0145
0146 std::vector<ME::Time> ME::timeDiff(Time t1, Time t2, short int& sign) {
0147   sign = 1;
0148   //  Time t1 = time_high( T1 );
0149   //  Time t2 = time_high( T2 );
0150   Time dt_s(0);
0151   if (t1 > t2) {
0152     dt_s = t1 - t2;
0153   } else {
0154     sign = -1;
0155     dt_s = t2 - t1;
0156   }
0157   Time dt_min = dt_s / 60;
0158   Time n_s = dt_s - dt_min * 60;
0159   Time dt_h = dt_min / 60;
0160   Time n_min = dt_min - dt_h * 60;
0161   Time dt_day = dt_h / 24;
0162   Time n_h = dt_h - dt_day * 24;
0163   Time n_day = dt_day;
0164
0165   std::vector<Time> vec_;
0166   vec_.push_back(n_day);
0167   vec_.push_back(n_h);
0168   vec_.push_back(n_min);
0169   vec_.push_back(n_s);
0170
0171   return vec_;
0172 }
0173
0174 float ME::timeDiff(Time t1, Time t0, int tunit) {
0175   float sign = 1.;
0176   Time dt(0);
0177   if (t1 > t0) {
0178     dt = t1 - t0;
0179   } else {
0180     sign = -1.;
0181     dt = t0 - t1;
0182   }
0183   float dt_f = ((float)dt) * sign;
0184   switch (tunit) {
0185     case iDay:
0186       return dt_f / 86400.;
0187       break;
0188     case iHour:
0189       return dt_f / 3600.;
0190       break;
0191     case iMinute:
0192       return dt_f / 60.;
0193       break;
0194     default:
0195       return dt_f;
0196   };
0197   return 0;
0198 }
0199
0200 ME::Time ME::time(float dt, Time t0, int tunit) {
0201   short int sign = 1;
0202   if (dt < 0)
0203     sign = -1;
0204   float t_ = sign * dt;
0205   switch (tunit) {
0206     case iDay:
0207       t_ *= 86400;
0208       break;
0209     case iHour:
0210       t_ *= 3600;
0211       break;
0212     case iMinute:
0213       t_ *= 60;
0214       break;
0215   };
0216   ME::Time it_ = static_cast<ME::Time>(t_);
0217   std::cout << "dt/it/t0/ " << dt << "/" << it_ << "/" << t0 << std::endl;
0218   if (sign == 1)
0219     return t0 + it_;
0220   else
0221     return t0 - it_;
0222 }
0223
0224 ME::Time ME::time_low(TimeStamp t) { return static_cast<Time>(kLowMask & t); }
0225
0226 ME::Time ME::time_high(TimeStamp t) { return static_cast<Time>(t >> 32); }
0227
0228 const TString ME::region[ME::iSizeE] = {"EE-", "EB-", "EB+", "EE+"};
0229
0230 int ME::ecalRegion(int ilmr) {
0231   assert(ilmr > 0 && ilmr <= 92);
0232   if (ilmr <= 36)
0233     return iEBM;
0234   ilmr -= 36;
0235   if (ilmr <= 36)
0236     return iEBP;
0237   ilmr -= 36;
0238   if (ilmr <= 10)
0239     return iEEP;
0240   return iEEM;
0241 }
0242
0243 int ME::lmr(int idcc, int side) {
0244   int ilmr = 0;
0245
0246   assert(side == 0 || side == 1);
0247
0248   if (idcc > 600)
0249     idcc -= 600;
0250   assert(idcc >= 1 && idcc <= 54);
0251   int ireg;
0252   if (idcc <= 9)
0253     ireg = iEEM;
0254   else {
0255     idcc -= 9;
0256     if (idcc <= 18)
0257       ireg = iEBM;
0258     else {
0259       idcc -= 18;
0260       if (idcc <= 18)
0261         ireg = iEBP;
0262       else {
0263         idcc -= 18;
0264         if (idcc <= 9)
0265           ireg = iEEP;
0266         else
0267           abort();
0268       }
0269     }
0270   }
0271   if (ireg == iEEM || ireg == iEEP) {
0272     if (side == 1 && idcc != 8) {
0273       return -1;
0274     }
0275     ilmr = idcc;
0276     if (idcc == 9)
0277       ilmr++;
0278     if (idcc == 8 && side == 1)
0279       ilmr++;
0280   } else if (ireg == iEBM || ireg == iEBP) {
0281     ilmr = 2 * (idcc - 1) + side + 1;
0282   } else
0283     abort();
0284
0285   if (ireg == iEBP)
0286     ilmr += 36;
0287   else if (ireg == iEEP)
0288     ilmr += 72;
0289   else if (ireg == iEEM)
0290     ilmr += 82;
0291
0292   return ilmr;
0293 }
0294
0295 std::pair<int, int> ME::dccAndSide(int ilmr) {
0296   int idcc = 0;
0297   int side = 0;
0298
0299   int ireg = ecalRegion(ilmr);
0300   if (ireg == iEEM)
0301     ilmr -= 82;
0302   else if (ireg == iEBP)
0303     ilmr -= 36;
0304   else if (ireg == iEEP)
0305     ilmr -= 72;
0306
0307   if (ireg == iEEM || ireg == iEEP) {
0308     assert(ilmr >= 1 && ilmr <= 10);
0309     side = 0;
0310     idcc = ilmr;
0311     if (ilmr >= 9)
0312       idcc--;
0313     if (ilmr == 9)
0314       side = 1;
0315   } else {
0316     assert(ilmr >= 1 && ilmr <= 36);
0317     idcc = (ilmr - 1) / 2 + 1;
0318     side = (ilmr - 1) % 2;
0319   }
0320
0321   if (ireg > iEEM)
0322     idcc += 9;
0323   if (ireg > iEBM)
0324     idcc += 18;
0325   if (ireg > iEBP)
0326     idcc += 18;
0327
0328   //  idcc += 600;
0329
0330   return std::pair<int, int>(idcc, side);
0331 }
0332
0333 void ME::regionAndSector(int ilmr, int& ireg, int& ism, int& idcc, int& side) {
0334   ireg = ecalRegion(ilmr);
0335
0336   std::pair<int, int> ipair_ = dccAndSide(ilmr);
0337   idcc = ipair_.first;
0338   side = ipair_.second;
0339
0340   ism = 0;
0341   if (ireg == iEEM || ireg == iEEP) {
0342     if (idcc > 600)
0343       idcc -= 600;  // also works with FEDids
0344     if (idcc >= 1 && idcc <= 9) {
0345       ism = 6 + idcc;
0346       if (ism > 9)
0347         ism -= 9;
0348       ism += 9;
0349     } else if (idcc >= 46 && idcc <= 54) {
0350       ism = idcc - 46 + 7;
0351       if (ism > 9)
0352         ism -= 9;
0353     } else
0354       abort();
0355   } else if (ireg == iEBM || ireg == iEBP) {
0356     if (idcc > 600)
0357       idcc -= 600;  // also works with FEDids
0358     assert(idcc >= 10 && idcc <= 45);
0359     ism = idcc - 9;
0360     if (ism > 18)
0361       ism -= 18;
0362     else
0363       ism += 18;
0364   } else
0365     abort();
0366 }
0367
0368 TString ME::smName(int ireg, int ism) {
0369   TString out;
0370   if (ireg == ME::iEEM || ireg == ME::iEEP) {
0371     assert(ism >= 1 && ism <= 18);
0372     out = "EE+";
0373     if (ireg == ME::iEEM)
0374       out = "EE-";
0375     if (ism > 9)
0376       ism -= 9;
0377     out += ism;
0378   } else if (ireg == ME::iEBM || ireg == ME::iEBP) {
0379     assert(ism >= 1 && ism <= 36);
0380     out = "EB+";
0381     if (ism > 18) {
0382       out = "EB-";
0383       ism -= 18;
0384     }
0385     out += ism;
0386   } else
0387     abort();
0388   return out;
0389 }
0390
0391 TString ME::smName(int ilmr) {
0392   TString out;
0393   int reg_(0);
0394   int sm_(0);
0395   int dcc_(0);
0396   int side_(0);
0397   ME::regionAndSector(ilmr, reg_, sm_, dcc_, side_);
0398   out = smName(reg_, sm_);
0399   return out;
0400 }
0401
0402 TString ME::smNameFromDcc(int idcc) {
0403   int ilmr = lmr(idcc, 0);
0404   return smName(ilmr);
0405 }
0406
0407 MEChannel* ME::lmrTree(int ilmr) { return regTree(ecalRegion(ilmr))->getDescendant(iLMRegion, ilmr); }
0408
0409 MEChannel* ME::regTree(int ireg) {
0410   assert(ireg >= iEEM && ireg <= iEEP);
0411   if (_trees[ireg] != nullptr)
0412     return _trees[ireg];
0413
0414   int iEcalRegion_ = ireg;
0415   int iSector_ = 0;
0416   int iLMRegion_ = 0;
0417   int iLMModule_ = 0;
0418   int iSuperCrystal_ = 0;
0419   int iCrystal_ = 0;
0420   MEChannel* leaf_(nullptr);
0421   MEChannel* tree_(nullptr);
0422
0423   if (iEcalRegion_ == iEBM || iEcalRegion_ == iEBP) {
0424     for (int isect = 1; isect <= 18; isect++) {
0425       iSector_ = isect;
0426       if (iEcalRegion_ == iEBM)
0427         iSector_ += 18;
0428       if (_trees[iEcalRegion_] == nullptr) {
0429         //        std::cout << "Building the tree of crystals -- "
0430         //         << ME::region[iEcalRegion_];
0431         _trees[iEcalRegion_] = new MEChannel(0, 0, iEcalRegion_, nullptr);
0432       }
0433       tree_ = _trees[iEcalRegion_];
0434       for (int iX = 0; iX < 17; iX++) {
0435         for (int iY = 0; iY < 4; iY++) {
0436           iSuperCrystal_ = MEEBGeom::tt_channel(iX, iY);
0437           iLMModule_ = MEEBGeom::lm_channel(iX, iY);
0438           for (int jx = 0; jx < 5; jx++) {
0439             for (int jy = 0; jy < 5; jy++) {
0440               int ix = 5 * iX + jx;
0441               int iy = 5 * iY + jy;
0442               if (useElectronicNumbering) {
0443                 iCrystal_ = MEEBGeom::electronic_channel(ix, iy);
0444               } else {
0445                 iCrystal_ = MEEBGeom::crystal_channel(ix, iy);
0446               }
0447               MEEBGeom::EtaPhiCoord globalCoord = MEEBGeom::globalCoord(iSector_, ix, iy);
0448               int ieta = globalCoord.first;
0449               int iphi = globalCoord.second;
0450               iLMRegion_ = MEEBGeom::lmr(ieta, iphi);
0451               leaf_ = tree_;
0452               leaf_ = leaf_->getDaughter(ieta, iphi, iSector_);
0453               leaf_ = leaf_->getDaughter(ieta, iphi, iLMRegion_);
0454               leaf_ = leaf_->getDaughter(ieta, iphi, iLMModule_);
0455               leaf_ = leaf_->getDaughter(ieta, iphi, iSuperCrystal_);
0456               leaf_->getDaughter(ieta, iphi, iCrystal_);
0457             }
0458           }
0459         }
0460       }
0461     }
0462   } else if (iEcalRegion_ == iEEM || iEcalRegion_ == iEEP) {
0463     int iz = 1;
0464     if (iEcalRegion_ == iEEM)
0465       iz = -1;
0466     if (_trees[iEcalRegion_] == nullptr) {
0467       //      std::cout << "Building the tree of crystals -- "
0468       //           << ME::region[iEcalRegion_];
0469       _trees[iEcalRegion_] = new MEChannel(0, 0, iEcalRegion_, nullptr);
0470     }
0471     tree_ = _trees[iEcalRegion_];
0472
0473     for (int ilmr = 72; ilmr <= 92; ilmr++)  // force the order of Monitoring Regions
0474     {
0475       if (ecalRegion(ilmr) != iEcalRegion_)
0476         continue;
0477       for (int ilmm = 1; ilmm <= 19; ilmm++)  // force the order of Monitoring Modules
0478       {
0479         for (int iXX = 1; iXX <= 10; iXX++) {
0480           for (int iside = 1; iside <= 2; iside++)  // symmetrize wrt y-axis
0481           {
0482             int iX = iXX;
0483             if (iside == 2)
0484               iX = 20 - iXX + 1;
0485             for (int iY = 1; iY <= 20; iY++) {
0486               //int iSector_   = MEEEGeom::sector( iX, iY );
0487               int iSector_ = MEEEGeom::sm(iX, iY, iz);
0488               if (iSector_ < 0)
0489                 continue;
0490               iLMRegion_ = MEEEGeom::lmr(iX, iY, iz);
0491               if (iLMRegion_ != ilmr)
0492                 continue;
0493               iLMModule_ = MEEEGeom::lmmod(iX, iY);
0494               if (iLMModule_ != ilmm)
0495                 continue;
0496               iSuperCrystal_ = MEEEGeom::sc(iX, iY);
0497
0498               for (int jxx = 1; jxx <= 5; jxx++) {
0499                 int jx = jxx;  // symmetrize...
0500                 if (iside == 2)
0501                   jx = 5 - jxx + 1;
0502
0503                 for (int jy = 1; jy <= 5; jy++) {
0504                   int ix = 5 * (iX - 1) + jx;
0505                   int iy = 5 * (iY - 1) + jy;
0506                   iCrystal_ = MEEEGeom::crystal(ix, iy);
0507                   if (iCrystal_ < 0)
0508                     continue;
0509                   leaf_ = tree_;
0510                   leaf_ = leaf_->getDaughter(ix, iy, iSector_);
0511                   leaf_ = leaf_->getDaughter(ix, iy, iLMRegion_);
0512                   leaf_ = leaf_->getDaughter(ix, iy, iLMModule_);
0513                   leaf_ = leaf_->getDaughter(ix, iy, iSuperCrystal_);
0514                   leaf_->getDaughter(ix, iy, iCrystal_);
0515                 }
0516               }
0517             }
0518           }
0519         }
0520       }
0521     }
0522   }
0523   //  std::cout << ".... done" << std::endl;
0524   return _trees[iEcalRegion_];
0525 }
0526
0527 bool ME::isBarrel(int ilmr) {
0528   int reg_ = ecalRegion(ilmr);
0529   if (reg_ == iEEM || reg_ == iEEP)
0530     return false;
0531   else if (reg_ == iEBM || reg_ == iEBP)
0532     return true;
0533   else
0534     abort();
0535   return true;
0536 }
0537
0538 std::pair<int, int> ME::memFromLmr(int ilmr) {
0539   if (isBarrel(ilmr))
0540     return MEEBGeom::memFromLmr(ilmr);
0541   else
0542     return MEEEGeom::memFromLmr(ilmr);
0543   return std::pair<int, int>();
0544 }
0545 std::vector<int> ME::apdRefChannels(int ilmmod, int ilmr) {
0546   if (isBarrel(ilmr))
0547     return MEEBGeom::apdRefChannels(ilmmod);
0548   else
0549     return MEEEGeom::apdRefChannels(ilmmod);
0550   return std::vector<int>();
0551 }
0552
0553 std::vector<int> ME::lmmodFromLmr(int ilmr) {
0554   if (isBarrel(ilmr))
0555     return MEEBGeom::lmmodFromLmr(ilmr);
0556   else
0557     return MEEEGeom::lmmodFromLmr(ilmr);
0558   return std::vector<int>();
0559 }
0560
0561 std::vector<int> ME::memFromDcc(int idcc) {
0562   std::vector<int> vec;
0563   for (int iside = 0; iside <= 1; iside++) {
0564     int ilmr = lmr(idcc, iside);
0565     if (ilmr < 0)
0566       continue;
0567     std::pair<int, int> mem_ = memFromLmr(ilmr);
0568     vec.push_back(mem_.first);
0569     vec.push_back(mem_.second);
0570   }
0571   return vec;
0572 }
0573
0574 std::vector<int> ME::lmmodFromDcc(int idcc) {
0575   std::vector<int> vec;
0576   for (int iside = 0; iside <= 1; iside++) {
0577     int ilmr = lmr(idcc, iside);
0578     if (ilmr < 0)
0579       continue;
0580     bool isBarrel_ = isBarrel(ilmr);
0581     std::vector<int> vec_ = lmmodFromLmr(ilmr);
0582     for (unsigned ii = 0; ii < vec_.size(); ii++) {
0583       int ilmmod_ = vec_[ii];
0584       if (!isBarrel_) {
0585         // special case for Julie
0586         if (ilmmod_ == 18 && iside == 1)
0587           ilmmod_ = 20;
0588         if (ilmmod_ == 19 && iside == 1)
0589           ilmmod_ = 21;
0590       }
0591       vec.push_back(ilmmod_);
0592     }
0593   }
0594   return vec;
0595 }
0596
0597 std::pair<int, int> ME::pn(int ilmr, int ilmmod, ME::PN ipn) {
0598   std::pair<int, int> pnpair_(0, 0);
0599   std::pair<int, int> mempair_ = memFromLmr(ilmr);
0600   if (isBarrel(ilmr)) {
0601     pnpair_ = MEEBGeom::pn(ilmmod);
0602   } else {
0603     int dee_ = MEEEGeom::dee(ilmr);
0604     pnpair_ = MEEEGeom::pn(dee_, ilmmod);
0605   }
0606   int mem_(0);
0607   int pn_(0);
0608   if (ipn == iPNA) {
0609     mem_ = mempair_.first;
0610     pn_ = pnpair_.first;
0611   } else {
0612     mem_ = mempair_.second;
0613     pn_ = pnpair_.second;
0614   }
0615   return std::pair<int, int>(mem_, pn_);
0616 }