Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:58 

0001 #include "RecoTBCalo/HcalTBObjectUnpacker/interface/HcalTBTDCUnpacker.h"
0002 #include "FWCore/Utilities/interface/Exception.h"
0003 #include <cmath>
0004 
0005 // Timing channels
0006 static const int lcTriggerTime = 1;
0007 static const int lcTTCL1ATime = 2;
0008 static const int lcBeamCoincidence = 3;
0009 static const int lcLaserFlash = 4;
0010 static const int lcQIEPhase = 5;
0011 
0012 static const int lcMuon1 = 90;
0013 static const int lcMuon2 = 91;
0014 static const int lcMuon3 = 92;
0015 static const int lcScint1 = 93;
0016 static const int lcScint2 = 94;
0017 static const int lcScint3 = 95;
0018 static const int lcScint4 = 96;
0019 static const int lcBeamHalo1 = 50;
0020 static const int lcBeamHalo2 = 51;
0021 static const int lcBeamHalo3 = 55;
0022 static const int lcBeamHalo4 = 63;
0023 static const int lcTOF1S = 129;
0024 static const int lcTOF1J = 130;
0025 static const int lcTOF2S = 131;
0026 static const int lcTOF2J = 132;
0027 
0028 namespace hcaltb {
0029 
0030   HcalTBTDCUnpacker::HcalTBTDCUnpacker(bool include_unmatched_hits)
0031       : includeUnmatchedHits_(include_unmatched_hits), dumpObs_(nullptr) {
0032     //  setupWC(); reads it from configuration file
0033     //  dumpObs_=fopen("dump_obs.csv","w");
0034   }
0035   void HcalTBTDCUnpacker::setCalib(const std::vector<std::vector<std::string> >& calibLines_) {
0036     for (int i = 0; i < 161; i++) {
0037       tdc_ped[i] = 0.;
0038       tdc_convers[i] = 1.;
0039     }
0040     for (unsigned int ii = 0; ii < calibLines_.size(); ii++) {
0041       //   TDC configuration
0042       if (calibLines_[ii][0] == "TDC") {
0043         if (calibLines_[ii].size() == 4) {
0044           int channel = atoi(calibLines_[ii][1].c_str());
0045           tdc_ped[channel] = atof(calibLines_[ii][2].c_str());
0046           tdc_convers[channel] = atof(calibLines_[ii][3].c_str());
0047           //        printf("Got TDC %i ped %f , conversion %f\n",channel, tdc_ped[channel],tdc_convers[channel]);
0048         } else {
0049           throw cms::Exception("Incomplete configuration")
0050               << "Wrong TDC configuration format : expected 3 parameters, got " << calibLines_[ii].size() - 1;
0051         }
0052       }  // End of the TDCs
0053 
0054       //   Wire chambers calibration
0055       if (calibLines_[ii][0] == "WC") {
0056         if (calibLines_[ii].size() == 6) {
0057           int plane = atoi(calibLines_[ii][1].c_str());
0058           wc_[plane].b0 = atof(calibLines_[ii][2].c_str());
0059           wc_[plane].b1 = atof(calibLines_[ii][3].c_str());
0060           wc_[plane].mean = atof(calibLines_[ii][4].c_str());
0061           wc_[plane].sigma = atof(calibLines_[ii][5].c_str());
0062           //         printf("Got WC plane %i b0 %f, b1 %f, mean %f, sigma %f\n",plane,
0063           //         wc_[plane].b0,wc_[plane].b1,wc_[plane].mean,wc_[plane].sigma);
0064         } else {
0065           throw cms::Exception("Incomplete configuration")
0066               << "Wrong Wire Chamber configuration format : expected 5 parameters, got " << calibLines_[ii].size() - 1;
0067         }
0068       }  // End of the Wire Chambers
0069 
0070     }  // End of the CalibLines
0071   }
0072 
0073   void HcalTBTDCUnpacker::unpack(const FEDRawData& raw, HcalTBEventPosition& pos, HcalTBTiming& timing) const {
0074     std::vector<Hit> hits;
0075 
0076     unpackHits(raw, hits, timing);
0077 
0078     reconstructWC(hits, pos);
0079     reconstructTiming(hits, timing);
0080   }
0081 
0082   struct ClassicTDCDataFormat {
0083     unsigned int cdfHeader0, cdfHeader1, cdfHeader2, cdfHeader3;
0084     unsigned int n_max_hits;  // maximum number of hits possible in the block
0085     unsigned int n_hits;
0086     unsigned int hits[2];
0087   };
0088 
0089   struct CombinedTDCQDCDataFormat {
0090     unsigned int cdfHeader0, cdfHeader1, cdfHeader2, cdfHeader3;
0091     unsigned int n_qdc_hits;  // Count of QDC channels
0092     unsigned int n_tdc_hits;  // upper/lower TDC counts
0093     unsigned short qdc_values[4];
0094   };
0095 
0096   //static const double CONVERSION_FACTOR=25.0/32.0;
0097 
0098   void HcalTBTDCUnpacker::unpackHits(const FEDRawData& raw, std::vector<Hit>& hits, HcalTBTiming& timing) const {
0099     const ClassicTDCDataFormat* tdc = (const ClassicTDCDataFormat*)raw.data();
0100 
0101     if (raw.size() < 3 * 8) {
0102       throw cms::Exception("Missing Data") << "No data in the TDC block";
0103     }
0104 
0105     const unsigned int* hitbase = nullptr;
0106     unsigned int totalhits = 0;
0107 
0108     // old TDC (767)
0109     if (tdc->n_max_hits != 192) {
0110       const CombinedTDCQDCDataFormat* qdctdc = (const CombinedTDCQDCDataFormat*)raw.data();
0111       hitbase = (const unsigned int*)(qdctdc);
0112       hitbase += 6;                             // header
0113       hitbase += qdctdc->n_qdc_hits / 2;        // two unsigned short per unsigned long
0114       totalhits = qdctdc->n_tdc_hits & 0xFFFF;  // mask off high bits
0115 
0116       //    for (unsigned int i=0; i<qdctdc->n_qdc_hits; i++)
0117       //      printf("QADC: %02d %d\n",i,qdctdc->qdc_values[i]&0xFFF);
0118 
0119     } else {
0120       hitbase = &(tdc->hits[0]);
0121       totalhits = tdc->n_hits;
0122     }
0123 
0124     for (unsigned int i = 0; i < totalhits; i++) {
0125       Hit h;
0126       h.channel = (hitbase[i] & 0x7FC00000) >> 22;  // hardcode channel assignment
0127       h.time = (hitbase[i] & 0xFFFFF) * tdc_convers[h.channel];
0128       hits.push_back(h);
0129       //        printf("V767: %d %f\n",h.channel,h.time);
0130     }
0131 
0132     // new TDC (V775)
0133     int v775[32];
0134     for (int i = 0; i < 32; i++)
0135       v775[i] = -1;
0136     if (tdc->n_max_hits != 192) {
0137       const CombinedTDCQDCDataFormat* qdctdc = (const CombinedTDCQDCDataFormat*)raw.data();
0138       hitbase = (const unsigned int*)(qdctdc);
0139       hitbase += 6;                                         // header
0140       hitbase += qdctdc->n_qdc_hits / 2;                    // two unsigned short per unsigned long
0141       hitbase += (qdctdc->n_tdc_hits & 0xFFFF);             // same length
0142       totalhits = (qdctdc->n_tdc_hits & 0xFFFF0000) >> 16;  // mask off high bits
0143       for (unsigned int i = 0; i < totalhits; i++) {
0144         Hit h;
0145         //      h.channel=129+i;
0146         h.channel = 129 + ((hitbase[i] & 0x3F0000) >> 16);
0147         h.time = (hitbase[i] & 0xFFF) * tdc_convers[h.channel];
0148         hits.push_back(h);
0149         if ((h.channel - 129) < 32)
0150           v775[(h.channel - 129)] = (hitbase[i] & 0xFFF);
0151         //      printf("V775: %d %f\n",h.channel,h.time);
0152       }
0153     }
0154     timing.setV775(v775);
0155   }
0156 
0157   void HcalTBTDCUnpacker::reconstructTiming(const std::vector<Hit>& hits, HcalTBTiming& timing) const {
0158     std::vector<Hit>::const_iterator j;
0159     double trigger_time = 0;
0160     double ttc_l1a_time = 0;
0161     double laser_flash = 0;
0162     double qie_phase = 0;
0163     double TOF1S_time = 0;
0164     double TOF1J_time = 0;
0165     double TOF2S_time = 0;
0166     double TOF2J_time = 0;
0167 
0168     std::vector<double> m1hits, m2hits, m3hits, s1hits, s2hits, s3hits, s4hits, bh1hits, bh2hits, bh3hits, bh4hits,
0169         beam_coinc;
0170 
0171     for (j = hits.begin(); j != hits.end(); j++) {
0172       switch (j->channel) {
0173         case lcTriggerTime:
0174           trigger_time = j->time - tdc_ped[lcTriggerTime];
0175           break;
0176         case lcTTCL1ATime:
0177           ttc_l1a_time = j->time - tdc_ped[lcTTCL1ATime];
0178           break;
0179         case lcBeamCoincidence:
0180           beam_coinc.push_back(j->time - tdc_ped[lcBeamCoincidence]);
0181           break;
0182         case lcLaserFlash:
0183           laser_flash = j->time - tdc_ped[lcLaserFlash];
0184           break;
0185         case lcQIEPhase:
0186           qie_phase = j->time - tdc_ped[lcQIEPhase];
0187           break;
0188         case lcMuon1:
0189           m1hits.push_back(j->time - tdc_ped[lcMuon1]);
0190           break;
0191         case lcMuon2:
0192           m2hits.push_back(j->time - tdc_ped[lcMuon2]);
0193           break;
0194         case lcMuon3:
0195           m3hits.push_back(j->time - tdc_ped[lcMuon3]);
0196           break;
0197         case lcScint1:
0198           s1hits.push_back(j->time - tdc_ped[lcScint1]);
0199           break;
0200         case lcScint2:
0201           s2hits.push_back(j->time - tdc_ped[lcScint2]);
0202           break;
0203         case lcScint3:
0204           s3hits.push_back(j->time - tdc_ped[lcScint3]);
0205           break;
0206         case lcScint4:
0207           s4hits.push_back(j->time - tdc_ped[lcScint4]);
0208           break;
0209         case lcTOF1S:
0210           TOF1S_time = j->time - tdc_ped[lcTOF1S];
0211           break;
0212         case lcTOF1J:
0213           TOF1J_time = j->time - tdc_ped[lcTOF1J];
0214           break;
0215         case lcTOF2S:
0216           TOF2S_time = j->time - tdc_ped[lcTOF2S];
0217           break;
0218         case lcTOF2J:
0219           TOF2J_time = j->time - tdc_ped[lcTOF2J];
0220           break;
0221         case lcBeamHalo1:
0222           bh1hits.push_back(j->time - tdc_ped[lcBeamHalo1]);
0223           break;
0224         case lcBeamHalo2:
0225           bh2hits.push_back(j->time - tdc_ped[lcBeamHalo2]);
0226           break;
0227         case lcBeamHalo3:
0228           bh3hits.push_back(j->time - tdc_ped[lcBeamHalo3]);
0229           break;
0230         case lcBeamHalo4:
0231           bh4hits.push_back(j->time - tdc_ped[lcBeamHalo4]);
0232           break;
0233         default:
0234           break;
0235       }
0236     }
0237 
0238     timing.setTimes(trigger_time, ttc_l1a_time, laser_flash, qie_phase, TOF1S_time, TOF1J_time, TOF2S_time, TOF2J_time);
0239     timing.setHits(
0240         m1hits, m2hits, m3hits, s1hits, s2hits, s3hits, s4hits, bh1hits, bh2hits, bh3hits, bh4hits, beam_coinc);
0241   }
0242 
0243   const int HcalTBTDCUnpacker::WC_CHANNELIDS[PLANECOUNT * 3] = {
0244       12,  13,  14,   // WCA LR plane
0245       10,  11,  14,   // WCA UD plane
0246       22,  23,  24,   // WCB LR plane
0247       20,  21,  24,   // WCB UD plane
0248       32,  33,  34,   // WCC LR plane
0249       30,  31,  34,   // WCC UD plane
0250       101, 102, 104,  // WCD LR plane
0251       107, 108, 110,  // WCD UD plane
0252       113, 114, 116,  // WCE LR plane
0253       97,  98,  99,   // WCE UD plane
0254       42,  43,  -1,   // WCF LR plane (was WC1)
0255       44,  60,  -1,   // WCF UD plane (was WC1)
0256       40,  41,  -1,   // WCG LR plane (was WC2)
0257       45,  61,  -1,   // WCG UD plane (was WC2)
0258       52,  53,  -1,   // WCH LR plane (was WC3)
0259       54,  62,  -1    // WCH UD plane (was WC3)
0260   };
0261 
0262   static const double TDC_OFFSET_CONSTANT = 12000;
0263   static const double N_SIGMA = 2.5;
0264 
0265   /* Obsolated - reads it from the configuration file
0266 void HcalTBTDCUnpacker::setupWC() {
0267 
0268   wc_[0].b0 = -0.0870056; wc_[0].b1 = -0.193263;  // WCA planes
0269   wc_[1].b0 = -0.0288171; wc_[1].b1 = -0.191231;
0270 
0271   wc_[2].b0 = -0.2214840; wc_[2].b1 = -0.191683;  // WCB planes
0272   wc_[3].b0 = -1.0847300; wc_[3].b1 = -0.187691;
0273 
0274   wc_[4].b0 = -0.1981440; wc_[4].b1 = -0.192760;  // WCC planes
0275   wc_[5].b0 =  0.4230690; wc_[5].b1 = -0.192278;
0276 
0277   wc_[6].b0 = -0.6039130; wc_[6].b1 = -0.185674;  // WCD planes
0278   wc_[7].b0 = -0.4366590; wc_[7].b1 = -0.184992;
0279 
0280   wc_[8].b0 =  1.7016400; wc_[8].b1 = -0.185575;  // WCE planes
0281   wc_[9].b0 = -0.2324480; wc_[9].b1 = -0.185367;
0282 
0283   wc_[0].mean=225.2; wc_[0].sigma=5.704; 
0284   wc_[1].mean=223.5; wc_[1].sigma=5.862; 
0285   wc_[2].mean=227.2; wc_[2].sigma=5.070; 
0286   wc_[3].mean=235.7; wc_[3].sigma=6.090; 
0287   wc_[4].mean=243.3; wc_[4].sigma=7.804; 
0288   wc_[5].mean=230.3; wc_[5].sigma=28.91; 
0289 
0290   wc_[6].mean=225.0; wc_[6].sigma=6.000; 
0291   wc_[7].mean=225.0; wc_[7].sigma=6.000; 
0292   wc_[8].mean=225.0; wc_[8].sigma=6.000; 
0293   wc_[9].mean=225.0; wc_[9].sigma=6.000; 
0294 }
0295 */
0296 
0297   void HcalTBTDCUnpacker::reconstructWC(const std::vector<Hit>& hits, HcalTBEventPosition& ep) const {
0298     // process all planes, looping over all hits...
0299     const int MAX_HITS = 100;
0300     float hits1[MAX_HITS], hits2[MAX_HITS], hitsA[MAX_HITS];
0301     int n1, n2, nA, chan1, chan2, chanA;
0302 
0303     std::vector<double> x;
0304 
0305     for (int plane = 0; plane < PLANECOUNT; plane++) {
0306       n1 = 0;
0307       n2 = 0;
0308       nA = 0;
0309 
0310       std::vector<double> plane_hits;
0311       double hit_time;
0312       hits1[0] = -1000;
0313       hits2[0] = -1000;
0314       hitsA[0] = -1000;
0315       hitsA[1] = -1000;
0316 
0317       chan1 = WC_CHANNELIDS[plane * 3];
0318       chan2 = WC_CHANNELIDS[plane * 3 + 1];
0319       chanA = WC_CHANNELIDS[plane * 3 + 2];
0320 
0321       for (std::vector<Hit>::const_iterator j = hits.begin(); j != hits.end(); j++) {
0322         if (j->channel == chan1 && n1 < MAX_HITS) {
0323           hits1[n1] = j->time - TDC_OFFSET_CONSTANT;
0324           n1++;
0325         }
0326         if (j->channel == chan2 && n2 < MAX_HITS) {
0327           hits2[n2] = j->time - TDC_OFFSET_CONSTANT;
0328           n2++;
0329         }
0330         if (j->channel == chanA && nA < MAX_HITS) {
0331           hitsA[nA] = j->time - TDC_OFFSET_CONSTANT;
0332           nA++;
0333         }
0334       }
0335 
0336       if (n1 != 0 && n2 != 0 && dumpObs_ != nullptr) {
0337         fprintf(dumpObs_, "%d,%f,%f,%f,%f\n", plane, hits1[0], hits2[0], hitsA[0], hitsA[1]);
0338         fflush(dumpObs_);
0339       }
0340 
0341       // anode-matched hits
0342       for (int ii = 0; ii < n1; ii++) {
0343         int jmin = -1, lmin = -1;
0344         float dsumMin = 99999;
0345         for (int jj = 0; jj < n2; jj++) {
0346           for (int ll = 0; ll < nA; ll++) {
0347             float dsum = fabs(wc_[plane].mean - hits1[ii] - hits2[jj] + 2.0 * hitsA[ll]);
0348             if (dsum < (N_SIGMA * wc_[plane].sigma) && dsum < dsumMin) {
0349               jmin = jj;
0350               lmin = ll;
0351               dsumMin = dsum;
0352             }
0353           }
0354         }
0355         if (jmin >= 0) {
0356           hit_time = wc_[plane].b0 + wc_[plane].b1 * (hits1[ii] - hits2[jmin]);
0357           if ((plane % 2) == 0) {
0358             plane_hits.push_back(-hit_time);
0359           } else {
0360             plane_hits.push_back(hit_time);
0361           }
0362           hits1[ii] = -99999;
0363           hits2[jmin] = -99999;
0364           hitsA[lmin] = 99999;
0365         }
0366       }
0367 
0368       if (includeUnmatchedHits_ || plane > 9)  // unmatched hits (all pairs get in here)
0369         for (int ii = 0; ii < n1; ii++) {
0370           if (hits1[ii] < -99990)
0371             continue;
0372           for (int jj = 0; jj < n2; jj++) {
0373             if (hits2[jj] < -99990)
0374               continue;
0375             hit_time = wc_[plane].b0 + wc_[plane].b1 * (hits1[ii] - hits2[jj]);
0376             if ((plane % 2) == 0) {
0377               plane_hits.push_back(-hit_time);
0378             } else {
0379               plane_hits.push_back(hit_time);
0380             }
0381           }
0382         }
0383 
0384       if ((plane % 2) == 0)
0385         x = plane_hits;
0386       else {
0387         char chamber = 'A' + plane / 2;
0388         ep.setChamberHits(chamber, x, plane_hits);
0389       }
0390     }
0391   }
0392 
0393 }  // namespace hcaltb
[ 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!