HcalTBTDCUnpacker.cc

Go to the documentation of this file.

#include "RecoTBCalo/HcalTBObjectUnpacker/interface/HcalTBTDCUnpacker.h"
#include "FWCore/Utilities/interface/Exception.h"
#include <cmath>

// Timing channels
static const int lcTriggerTime = 1;
static const int lcTTCL1ATime = 2;
static const int lcBeamCoincidence = 3;
static const int lcLaserFlash = 4;
static const int lcQIEPhase = 5;

static const int lcMuon1 = 90;
static const int lcMuon2 = 91;
static const int lcMuon3 = 92;
static const int lcScint1 = 93;
static const int lcScint2 = 94;
static const int lcScint3 = 95;
static const int lcScint4 = 96;
static const int lcBeamHalo1 = 50;
static const int lcBeamHalo2 = 51;
static const int lcBeamHalo3 = 55;
static const int lcBeamHalo4 = 63;
static const int lcTOF1S = 129;
static const int lcTOF1J = 130;
static const int lcTOF2S = 131;
static const int lcTOF2J = 132;

namespace hcaltb {

  HcalTBTDCUnpacker::HcalTBTDCUnpacker(bool include_unmatched_hits)
      : includeUnmatchedHits_(include_unmatched_hits), dumpObs_(nullptr) {
    //  setupWC(); reads it from configuration file
    //  dumpObs_=fopen("dump_obs.csv","w");
  }
  void HcalTBTDCUnpacker::setCalib(const std::vector<std::vector<std::string> >& calibLines_) {
    for (int i = 0; i < 161; i++) {
      tdc_ped[i] = 0.;
      tdc_convers[i] = 1.;
    }
    for (unsigned int ii = 0; ii < calibLines_.size(); ii++) {
      //   TDC configuration
      if (calibLines_[ii][0] == "TDC") {
        if (calibLines_[ii].size() == 4) {
          int channel = atoi(calibLines_[ii][1].c_str());
          tdc_ped[channel] = atof(calibLines_[ii][2].c_str());
          tdc_convers[channel] = atof(calibLines_[ii][3].c_str());
          //        printf("Got TDC %i ped %f , conversion %f\n",channel, tdc_ped[channel],tdc_convers[channel]);
        } else {
          throw cms::Exception("Incomplete configuration")
              << "Wrong TDC configuration format : expected 3 parameters, got " << calibLines_[ii].size() - 1;
        }
      }  // End of the TDCs

      //   Wire chambers calibration
      if (calibLines_[ii][0] == "WC") {
        if (calibLines_[ii].size() == 6) {
          int plane = atoi(calibLines_[ii][1].c_str());
          wc_[plane].b0 = atof(calibLines_[ii][2].c_str());
          wc_[plane].b1 = atof(calibLines_[ii][3].c_str());
          wc_[plane].mean = atof(calibLines_[ii][4].c_str());
          wc_[plane].sigma = atof(calibLines_[ii][5].c_str());
          //         printf("Got WC plane %i b0 %f, b1 %f, mean %f, sigma %f\n",plane,
          //         wc_[plane].b0,wc_[plane].b1,wc_[plane].mean,wc_[plane].sigma);
        } else {
          throw cms::Exception("Incomplete configuration")
              << "Wrong Wire Chamber configuration format : expected 5 parameters, got " << calibLines_[ii].size() - 1;
        }
      }  // End of the Wire Chambers

    }  // End of the CalibLines
  }

  void HcalTBTDCUnpacker::unpack(const FEDRawData& raw, HcalTBEventPosition& pos, HcalTBTiming& timing) const {
    std::vector<Hit> hits;

    unpackHits(raw, hits, timing);

    reconstructWC(hits, pos);
    reconstructTiming(hits, timing);
  }

  struct ClassicTDCDataFormat {
    unsigned int cdfHeader0, cdfHeader1, cdfHeader2, cdfHeader3;
    unsigned int n_max_hits;  // maximum number of hits possible in the block
    unsigned int n_hits;
    unsigned int hits[2];
  };

  struct CombinedTDCQDCDataFormat {
    unsigned int cdfHeader0, cdfHeader1, cdfHeader2, cdfHeader3;
    unsigned int n_qdc_hits;  // Count of QDC channels
    unsigned int n_tdc_hits;  // upper/lower TDC counts
    unsigned short qdc_values[4];
  };

  //static const double CONVERSION_FACTOR=25.0/32.0;

  void HcalTBTDCUnpacker::unpackHits(const FEDRawData& raw, std::vector<Hit>& hits, HcalTBTiming& timing) const {
    const ClassicTDCDataFormat* tdc = (const ClassicTDCDataFormat*)raw.data();

    if (raw.size() < 3 * 8) {
      throw cms::Exception("Missing Data") << "No data in the TDC block";
    }

    const unsigned int* hitbase = nullptr;
    unsigned int totalhits = 0;

    // old TDC (767)
    if (tdc->n_max_hits != 192) {
      const CombinedTDCQDCDataFormat* qdctdc = (const CombinedTDCQDCDataFormat*)raw.data();
      hitbase = (const unsigned int*)(qdctdc);
      hitbase += 6;                             // header
      hitbase += qdctdc->n_qdc_hits / 2;        // two unsigned short per unsigned long
      totalhits = qdctdc->n_tdc_hits & 0xFFFF;  // mask off high bits

      //    for (unsigned int i=0; i<qdctdc->n_qdc_hits; i++)
      //      printf("QADC: %02d %d\n",i,qdctdc->qdc_values[i]&0xFFF);

    } else {
      hitbase = &(tdc->hits[0]);
      totalhits = tdc->n_hits;
    }

    for (unsigned int i = 0; i < totalhits; i++) {
      Hit h;
      h.channel = (hitbase[i] & 0x7FC00000) >> 22;  // hardcode channel assignment
      h.time = (hitbase[i] & 0xFFFFF) * tdc_convers[h.channel];
      hits.push_back(h);
      //        printf("V767: %d %f\n",h.channel,h.time);
    }

    // new TDC (V775)
    int v775[32];
    for (int i = 0; i < 32; i++)
      v775[i] = -1;
    if (tdc->n_max_hits != 192) {
      const CombinedTDCQDCDataFormat* qdctdc = (const CombinedTDCQDCDataFormat*)raw.data();
      hitbase = (const unsigned int*)(qdctdc);
      hitbase += 6;                                         // header
      hitbase += qdctdc->n_qdc_hits / 2;                    // two unsigned short per unsigned long
      hitbase += (qdctdc->n_tdc_hits & 0xFFFF);             // same length
      totalhits = (qdctdc->n_tdc_hits & 0xFFFF0000) >> 16;  // mask off high bits
      for (unsigned int i = 0; i < totalhits; i++) {
        Hit h;
        //      h.channel=129+i;
        h.channel = 129 + ((hitbase[i] & 0x3F0000) >> 16);
        h.time = (hitbase[i] & 0xFFF) * tdc_convers[h.channel];
        hits.push_back(h);
        if ((h.channel - 129) < 32)
          v775[(h.channel - 129)] = (hitbase[i] & 0xFFF);
        //      printf("V775: %d %f\n",h.channel,h.time);
      }
    }
    timing.setV775(v775);
  }

  void HcalTBTDCUnpacker::reconstructTiming(const std::vector<Hit>& hits, HcalTBTiming& timing) const {
    std::vector<Hit>::const_iterator j;
    double trigger_time = 0;
    double ttc_l1a_time = 0;
    double laser_flash = 0;
    double qie_phase = 0;
    double TOF1S_time = 0;
    double TOF1J_time = 0;
    double TOF2S_time = 0;
    double TOF2J_time = 0;

    std::vector<double> m1hits, m2hits, m3hits, s1hits, s2hits, s3hits, s4hits, bh1hits, bh2hits, bh3hits, bh4hits,
        beam_coinc;

    for (j = hits.begin(); j != hits.end(); j++) {
      switch (j->channel) {
        case lcTriggerTime:
          trigger_time = j->time - tdc_ped[lcTriggerTime];
          break;
        case lcTTCL1ATime:
          ttc_l1a_time = j->time - tdc_ped[lcTTCL1ATime];
          break;
        case lcBeamCoincidence:
          beam_coinc.push_back(j->time - tdc_ped[lcBeamCoincidence]);
          break;
        case lcLaserFlash:
          laser_flash = j->time - tdc_ped[lcLaserFlash];
          break;
        case lcQIEPhase:
          qie_phase = j->time - tdc_ped[lcQIEPhase];
          break;
        case lcMuon1:
          m1hits.push_back(j->time - tdc_ped[lcMuon1]);
          break;
        case lcMuon2:
          m2hits.push_back(j->time - tdc_ped[lcMuon2]);
          break;
        case lcMuon3:
          m3hits.push_back(j->time - tdc_ped[lcMuon3]);
          break;
        case lcScint1:
          s1hits.push_back(j->time - tdc_ped[lcScint1]);
          break;
        case lcScint2:
          s2hits.push_back(j->time - tdc_ped[lcScint2]);
          break;
        case lcScint3:
          s3hits.push_back(j->time - tdc_ped[lcScint3]);
          break;
        case lcScint4:
          s4hits.push_back(j->time - tdc_ped[lcScint4]);
          break;
        case lcTOF1S:
          TOF1S_time = j->time - tdc_ped[lcTOF1S];
          break;
        case lcTOF1J:
          TOF1J_time = j->time - tdc_ped[lcTOF1J];
          break;
        case lcTOF2S:
          TOF2S_time = j->time - tdc_ped[lcTOF2S];
          break;
        case lcTOF2J:
          TOF2J_time = j->time - tdc_ped[lcTOF2J];
          break;
        case lcBeamHalo1:
          bh1hits.push_back(j->time - tdc_ped[lcBeamHalo1]);
          break;
        case lcBeamHalo2:
          bh2hits.push_back(j->time - tdc_ped[lcBeamHalo2]);
          break;
        case lcBeamHalo3:
          bh3hits.push_back(j->time - tdc_ped[lcBeamHalo3]);
          break;
        case lcBeamHalo4:
          bh4hits.push_back(j->time - tdc_ped[lcBeamHalo4]);
          break;
        default:
          break;
      }
    }

    timing.setTimes(trigger_time, ttc_l1a_time, laser_flash, qie_phase, TOF1S_time, TOF1J_time, TOF2S_time, TOF2J_time);
    timing.setHits(
        m1hits, m2hits, m3hits, s1hits, s2hits, s3hits, s4hits, bh1hits, bh2hits, bh3hits, bh4hits, beam_coinc);
  }

  const int HcalTBTDCUnpacker::WC_CHANNELIDS[PLANECOUNT * 3] = {
      12,  13,  14,   // WCA LR plane
      10,  11,  14,   // WCA UD plane
      22,  23,  24,   // WCB LR plane
      20,  21,  24,   // WCB UD plane
      32,  33,  34,   // WCC LR plane
      30,  31,  34,   // WCC UD plane
      101, 102, 104,  // WCD LR plane
      107, 108, 110,  // WCD UD plane
      113, 114, 116,  // WCE LR plane
      97,  98,  99,   // WCE UD plane
      42,  43,  -1,   // WCF LR plane (was WC1)
      44,  60,  -1,   // WCF UD plane (was WC1)
      40,  41,  -1,   // WCG LR plane (was WC2)
      45,  61,  -1,   // WCG UD plane (was WC2)
      52,  53,  -1,   // WCH LR plane (was WC3)
      54,  62,  -1    // WCH UD plane (was WC3)
  };

  static const double TDC_OFFSET_CONSTANT = 12000;
  static const double N_SIGMA = 2.5;

  /* Obsolated - reads it from the configuration file
void HcalTBTDCUnpacker::setupWC() {

  wc_[0].b0 = -0.0870056; wc_[0].b1 = -0.193263;  // WCA planes
  wc_[1].b0 = -0.0288171; wc_[1].b1 = -0.191231;

  wc_[2].b0 = -0.2214840; wc_[2].b1 = -0.191683;  // WCB planes
  wc_[3].b0 = -1.0847300; wc_[3].b1 = -0.187691;

  wc_[4].b0 = -0.1981440; wc_[4].b1 = -0.192760;  // WCC planes
  wc_[5].b0 =  0.4230690; wc_[5].b1 = -0.192278;

  wc_[6].b0 = -0.6039130; wc_[6].b1 = -0.185674;  // WCD planes
  wc_[7].b0 = -0.4366590; wc_[7].b1 = -0.184992;

  wc_[8].b0 =  1.7016400; wc_[8].b1 = -0.185575;  // WCE planes
  wc_[9].b0 = -0.2324480; wc_[9].b1 = -0.185367;

  wc_[0].mean=225.2; wc_[0].sigma=5.704; 
  wc_[1].mean=223.5; wc_[1].sigma=5.862; 
  wc_[2].mean=227.2; wc_[2].sigma=5.070; 
  wc_[3].mean=235.7; wc_[3].sigma=6.090; 
  wc_[4].mean=243.3; wc_[4].sigma=7.804; 
  wc_[5].mean=230.3; wc_[5].sigma=28.91; 

  wc_[6].mean=225.0; wc_[6].sigma=6.000; 
  wc_[7].mean=225.0; wc_[7].sigma=6.000; 
  wc_[8].mean=225.0; wc_[8].sigma=6.000; 
  wc_[9].mean=225.0; wc_[9].sigma=6.000; 
}
*/

  void HcalTBTDCUnpacker::reconstructWC(const std::vector<Hit>& hits, HcalTBEventPosition& ep) const {
    // process all planes, looping over all hits...
    const int MAX_HITS = 100;
    float hits1[MAX_HITS], hits2[MAX_HITS], hitsA[MAX_HITS];
    int n1, n2, nA, chan1, chan2, chanA;

    std::vector<double> x;

    for (int plane = 0; plane < PLANECOUNT; plane++) {
      n1 = 0;
      n2 = 0;
      nA = 0;

      std::vector<double> plane_hits;
      double hit_time;
      hits1[0] = -1000;
      hits2[0] = -1000;
      hitsA[0] = -1000;
      hitsA[1] = -1000;

      chan1 = WC_CHANNELIDS[plane * 3];
      chan2 = WC_CHANNELIDS[plane * 3 + 1];
      chanA = WC_CHANNELIDS[plane * 3 + 2];

      for (std::vector<Hit>::const_iterator j = hits.begin(); j != hits.end(); j++) {
        if (j->channel == chan1 && n1 < MAX_HITS) {
          hits1[n1] = j->time - TDC_OFFSET_CONSTANT;
          n1++;
        }
        if (j->channel == chan2 && n2 < MAX_HITS) {
          hits2[n2] = j->time - TDC_OFFSET_CONSTANT;
          n2++;
        }
        if (j->channel == chanA && nA < MAX_HITS) {
          hitsA[nA] = j->time - TDC_OFFSET_CONSTANT;
          nA++;
        }
      }

      if (n1 != 0 && n2 != 0 && dumpObs_ != nullptr) {
        fprintf(dumpObs_, "%d,%f,%f,%f,%f\n", plane, hits1[0], hits2[0], hitsA[0], hitsA[1]);
        fflush(dumpObs_);
      }

      // anode-matched hits
      for (int ii = 0; ii < n1; ii++) {
        int jmin = -1, lmin = -1;
        float dsumMin = 99999;
        for (int jj = 0; jj < n2; jj++) {
          for (int ll = 0; ll < nA; ll++) {
            float dsum = fabs(wc_[plane].mean - hits1[ii] - hits2[jj] + 2.0 * hitsA[ll]);
            if (dsum < (N_SIGMA * wc_[plane].sigma) && dsum < dsumMin) {
              jmin = jj;
              lmin = ll;
              dsumMin = dsum;
            }
          }
        }
        if (jmin >= 0) {
          hit_time = wc_[plane].b0 + wc_[plane].b1 * (hits1[ii] - hits2[jmin]);
          if ((plane % 2) == 0) {
            plane_hits.push_back(-hit_time);
          } else {
            plane_hits.push_back(hit_time);
          }
          hits1[ii] = -99999;
          hits2[jmin] = -99999;
          hitsA[lmin] = 99999;
        }
      }

      if (includeUnmatchedHits_ || plane > 9)  // unmatched hits (all pairs get in here)
        for (int ii = 0; ii < n1; ii++) {
          if (hits1[ii] < -99990)
            continue;
          for (int jj = 0; jj < n2; jj++) {
            if (hits2[jj] < -99990)
              continue;
            hit_time = wc_[plane].b0 + wc_[plane].b1 * (hits1[ii] - hits2[jj]);
            if ((plane % 2) == 0) {
              plane_hits.push_back(-hit_time);
            } else {
              plane_hits.push_back(hit_time);
            }
          }
        }

      if ((plane % 2) == 0)
        x = plane_hits;
      else {
        char chamber = 'A' + plane / 2;
        ep.setChamberHits(chamber, x, plane_hits);
      }
    }
  }

}  // namespace hcaltb