CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_2_7_hltpatch1/src/RecoTBCalo/HcalTBObjectUnpacker/src/HcalTBTDCUnpacker.cc

Go to the documentation of this file.
00001 #include "RecoTBCalo/HcalTBObjectUnpacker/interface/HcalTBTDCUnpacker.h"
00002 #include "FWCore/Utilities/interface/Exception.h"
00003 #include <math.h>
00004 
00005 // Timing channels
00006 static const int lcTriggerTime     = 1;
00007 static const int lcTTCL1ATime      = 2;
00008 static const int lcBeamCoincidence = 3;
00009 static const int lcLaserFlash      = 4;
00010 static const int lcQIEPhase        = 5;
00011 
00012 static const int lcMuon1           = 90;
00013 static const int lcMuon2           = 91;
00014 static const int lcMuon3           = 92;
00015 static const int lcScint1          = 93;
00016 static const int lcScint2          = 94;
00017 static const int lcScint3          = 95;
00018 static const int lcScint4          = 96;
00019 static const int lcBeamHalo1       = 50;
00020 static const int lcBeamHalo2       = 51;
00021 static const int lcBeamHalo3       = 55;
00022 static const int lcBeamHalo4       = 63;
00023 static const int lcTOF1S           = 129;
00024 static const int lcTOF1J           = 130;
00025 static const int lcTOF2S           = 131;
00026 static const int lcTOF2J           = 132;
00027 
00028 namespace hcaltb {
00029 
00030 HcalTBTDCUnpacker::HcalTBTDCUnpacker(bool include_unmatched_hits) :
00031   includeUnmatchedHits_(include_unmatched_hits) {
00032 //  setupWC(); reads it from configuration file
00033 }
00034 void HcalTBTDCUnpacker::setCalib(const std::vector<std::vector<std::string> >& calibLines_) {
00035         for(int i=0;i<161;i++)
00036          {
00037           tdc_ped[i]=0.;tdc_convers[i]=1.;
00038          }
00039         for(unsigned int ii=0;ii<calibLines_.size();ii++)
00040          {
00041 //   TDC configuration
00042           if(calibLines_[ii][0]=="TDC")
00043                 {
00044                 if(calibLines_[ii].size()==4)
00045                   {
00046                   int channel=atoi(calibLines_[ii][1].c_str());
00047                   tdc_ped[channel]=atof(calibLines_[ii][2].c_str());
00048                   tdc_convers[channel]=atof(calibLines_[ii][3].c_str());
00049         //        printf("Got TDC %i ped %f , conversion %f\n",channel, tdc_ped[channel],tdc_convers[channel]);
00050                   }
00051                  else
00052                   {
00053                throw cms::Exception("Incomplete configuration") << 
00054                 "Wrong TDC configuration format : expected 3 parameters, got "<<calibLines_[ii].size()-1;
00055                   }
00056                 } // End of the TDCs
00057 
00058 //   Wire chambers calibration
00059           if(calibLines_[ii][0]=="WC")
00060                 {
00061                 if(calibLines_[ii].size()==6)
00062                   {
00063                   int plane=atoi(calibLines_[ii][1].c_str());
00064                   wc_[plane].b0=atof(calibLines_[ii][2].c_str());
00065                   wc_[plane].b1=atof(calibLines_[ii][3].c_str());
00066                   wc_[plane].mean=atof(calibLines_[ii][4].c_str());
00067                   wc_[plane].sigma=atof(calibLines_[ii][5].c_str());
00068        //         printf("Got WC plane %i b0 %f, b1 %f, mean %f, sigma %f\n",plane, 
00069        //                wc_[plane].b0,wc_[plane].b1,wc_[plane].mean,wc_[plane].sigma);
00070                   }
00071                  else
00072                   {
00073                throw cms::Exception("Incomplete configuration") << 
00074                 "Wrong Wire Chamber configuration format : expected 5 parameters, got "<<calibLines_[ii].size()-1;
00075                   }
00076                 } // End of the Wire Chambers
00077 
00078          } // End of the CalibLines
00079         }
00080 
00081   void HcalTBTDCUnpacker::unpack(const FEDRawData& raw,
00082                                HcalTBEventPosition& pos,
00083                                HcalTBTiming& timing) const {
00084     std::vector<Hit> hits;
00085 
00086     unpackHits(raw, hits, timing);
00087     
00088     reconstructWC(hits, pos);
00089     reconstructTiming(hits, timing);
00090     
00091   }
00092 
00093   struct ClassicTDCDataFormat {
00094     unsigned int cdfHeader0,cdfHeader1,cdfHeader2,cdfHeader3;
00095     unsigned int n_max_hits; // maximum number of hits possible in the block
00096     unsigned int n_hits;
00097     unsigned int hits[2];
00098   };
00099 
00100   struct CombinedTDCQDCDataFormat {
00101     unsigned int cdfHeader0,cdfHeader1,cdfHeader2,cdfHeader3;
00102     unsigned int n_qdc_hits; // Count of QDC channels
00103     unsigned int n_tdc_hits; // upper/lower TDC counts    
00104     unsigned short qdc_values[4];
00105   };
00106 
00107 //static const double CONVERSION_FACTOR=25.0/32.0;
00108 
00109 void HcalTBTDCUnpacker::unpackHits(const FEDRawData& raw,
00110                                    std::vector<Hit>& hits,HcalTBTiming& timing) const {
00111   const ClassicTDCDataFormat* tdc=(const ClassicTDCDataFormat*)raw.data();
00112 
00113   if (raw.size()<3*8) {
00114     throw cms::Exception("Missing Data") << "No data in the TDC block";
00115   }
00116 
00117   const unsigned int* hitbase=0;
00118   unsigned int totalhits=0;
00119 
00120   // old TDC (767)
00121   if (tdc->n_max_hits!=192) {
00122     const CombinedTDCQDCDataFormat* qdctdc=(const CombinedTDCQDCDataFormat*)raw.data();
00123     hitbase=(unsigned int*)(qdctdc);
00124     hitbase+=6; // header
00125     hitbase+=qdctdc->n_qdc_hits/2; // two unsigned short per unsigned long
00126     totalhits=qdctdc->n_tdc_hits&0xFFFF; // mask off high bits
00127 
00128     //    for (unsigned int i=0; i<qdctdc->n_qdc_hits; i++)
00129     //      printf("QADC: %02d %d\n",i,qdctdc->qdc_values[i]&0xFFF);
00130 
00131   } else {
00132     hitbase=&(tdc->hits[0]);
00133     totalhits=tdc->n_hits;
00134   }
00135 
00136   for (unsigned int i=0; i<totalhits; i++) {
00137     Hit h;    
00138     h.channel=(hitbase[i]&0x7FC00000)>>22; // hardcode channel assignment
00139     h.time=(hitbase[i]&0xFFFFF)*tdc_convers[h.channel]; 
00140     hits.push_back(h);
00141     //        printf("V767: %d %f\n",h.channel,h.time);
00142   }
00143 
00144   // new TDC (V775)
00145   int v775[32];
00146   for (int i=0;i<32;i++) v775[i]=-1;
00147   if (tdc->n_max_hits!=192) {
00148     const CombinedTDCQDCDataFormat* qdctdc=(const CombinedTDCQDCDataFormat*)raw.data();
00149     hitbase=(unsigned int*)(qdctdc);
00150     hitbase+=6; // header
00151     hitbase+=qdctdc->n_qdc_hits/2; // two unsigned short per unsigned long
00152     hitbase+=(qdctdc->n_tdc_hits&0xFFFF); // same length
00153     totalhits=(qdctdc->n_tdc_hits&0xFFFF0000)>>16; // mask off high bits    
00154     for (unsigned int i=0; i<totalhits; i++) {
00155       Hit h;    
00156 //      h.channel=129+i;
00157       h.channel=129+((hitbase[i]&0x3F0000)>>16);
00158       h.time=(hitbase[i]&0xFFF)*tdc_convers[h.channel] ;
00159       hits.push_back(h);
00160       if ( (h.channel-129)<32 ) 
00161          v775[(h.channel-129)] = (hitbase[i]&0xFFF);
00162       //      printf("V775: %d %f\n",h.channel,h.time);
00163     }
00164   }
00165   timing.setV775(v775);
00166 }
00167 
00168 void HcalTBTDCUnpacker::reconstructTiming(const std::vector<Hit>& hits,
00169                                           HcalTBTiming& timing) const {
00170   std::vector<Hit>::const_iterator j;
00171   double trigger_time=0;
00172   double ttc_l1a_time=0;
00173   double laser_flash=0;
00174   double qie_phase=0;
00175   double TOF1S_time=0;
00176   double TOF1J_time=0;
00177   double TOF2S_time=0;
00178   double TOF2J_time=0;
00179   
00180   std::vector<double> m1hits, m2hits, m3hits, s1hits, s2hits, s3hits, s4hits,
00181                       bh1hits, bh2hits, bh3hits, bh4hits,beam_coinc;
00182 
00183   for (j=hits.begin(); j!=hits.end(); j++) {
00184     switch (j->channel) {
00185     case lcTriggerTime:     trigger_time   = j->time-tdc_ped[lcTriggerTime];  break;
00186     case lcTTCL1ATime:      ttc_l1a_time   = j->time-tdc_ped[lcTTCL1ATime];  break;
00187     case lcBeamCoincidence: beam_coinc.push_back(j->time-tdc_ped[lcBeamCoincidence]);  break;
00188     case lcLaserFlash:      laser_flash    = j->time-tdc_ped[lcLaserFlash];  break;
00189     case lcQIEPhase:        qie_phase      = j->time-tdc_ped[lcQIEPhase];  break;
00190     case lcMuon1:           m1hits.push_back(j->time-tdc_ped[lcMuon1]); break;
00191     case lcMuon2:           m2hits.push_back(j->time-tdc_ped[lcMuon2]); break;
00192     case lcMuon3:           m3hits.push_back(j->time-tdc_ped[lcMuon3]); break;
00193     case lcScint1:          s1hits.push_back(j->time-tdc_ped[lcScint1]); break;
00194     case lcScint2:          s2hits.push_back(j->time-tdc_ped[lcScint2]); break;
00195     case lcScint3:          s3hits.push_back(j->time-tdc_ped[lcScint3]); break;
00196     case lcScint4:          s4hits.push_back(j->time-tdc_ped[lcScint4]); break;
00197     case lcTOF1S:           TOF1S_time   = j->time-tdc_ped[lcTOF1S];  break;
00198     case lcTOF1J:           TOF1J_time   = j->time-tdc_ped[lcTOF1J];  break;
00199     case lcTOF2S:           TOF2S_time   = j->time-tdc_ped[lcTOF2S];  break;
00200     case lcTOF2J:           TOF2J_time   = j->time-tdc_ped[lcTOF2J];  break;
00201     case lcBeamHalo1:       bh1hits.push_back(j->time-tdc_ped[lcBeamHalo1]); break;
00202     case lcBeamHalo2:       bh2hits.push_back(j->time-tdc_ped[lcBeamHalo2]); break;
00203     case lcBeamHalo3:       bh3hits.push_back(j->time-tdc_ped[lcBeamHalo3]); break;
00204     case lcBeamHalo4:       bh4hits.push_back(j->time-tdc_ped[lcBeamHalo4]); break;
00205     default: break;
00206     }
00207   }
00208 
00209   timing.setTimes(trigger_time,ttc_l1a_time,laser_flash,qie_phase,TOF1S_time,TOF1J_time,TOF2S_time,TOF2J_time);
00210   timing.setHits (m1hits,m2hits,m3hits,s1hits,s2hits,s3hits,s4hits,bh1hits,bh2hits,bh3hits,bh4hits,beam_coinc);
00211 
00212 }
00213 
00214 const int HcalTBTDCUnpacker::WC_CHANNELIDS[PLANECOUNT*3] = { 
00215                                                      12, 13, 14, // WCA LR plane 
00216                                                      10, 11, 15, // WCA UD plane
00217                                                      22, 23, 24, // WCB LR plane
00218                                                      20, 21, 25, // WCB UD plane
00219                                                      32, 33, 34, // WCC LR plane
00220                                                      30, 31, 35, // WCC UD plane
00221                                                      101, 102, 104, // WCD LR plane
00222                                                      107, 108, 110, // WCD UD plane
00223                                                      113, 114, 116, // WCE LR plane
00224                                                      97, 98, 99, // WCE UD plane 
00225                                                     42, 43, -1, // WCF LR plane (was WC1)
00226                                                     44, 60, -1, // WCF UD plane (was WC1)
00227                                                     40, 41, -1, // WCG LR plane (was WC2)
00228                                                     45, 61, -1, // WCG UD plane (was WC2)
00229                                                     52, 53, -1, // WCH LR plane (was WC3)
00230                                                     54, 62, -1  // WCH UD plane (was WC3)
00231 };
00232 
00233 static const double TDC_OFFSET_CONSTANT = 12000;
00234 static const double N_SIGMA = 2.5;
00235 
00236 /* Obsolated - reads it from the configuration file
00237 void HcalTBTDCUnpacker::setupWC() {
00238 
00239   wc_[0].b0 = -0.0870056; wc_[0].b1 = -0.193263;  // WCA planes
00240   wc_[1].b0 = -0.0288171; wc_[1].b1 = -0.191231;
00241 
00242   wc_[2].b0 = -0.2214840; wc_[2].b1 = -0.191683;  // WCB planes
00243   wc_[3].b0 = -1.0847300; wc_[3].b1 = -0.187691;
00244 
00245   wc_[4].b0 = -0.1981440; wc_[4].b1 = -0.192760;  // WCC planes
00246   wc_[5].b0 =  0.4230690; wc_[5].b1 = -0.192278;
00247 
00248   wc_[6].b0 = -0.6039130; wc_[6].b1 = -0.185674;  // WCD planes
00249   wc_[7].b0 = -0.4366590; wc_[7].b1 = -0.184992;
00250 
00251   wc_[8].b0 =  1.7016400; wc_[8].b1 = -0.185575;  // WCE planes
00252   wc_[9].b0 = -0.2324480; wc_[9].b1 = -0.185367;
00253 
00254   wc_[0].mean=225.2; wc_[0].sigma=5.704; 
00255   wc_[1].mean=223.5; wc_[1].sigma=5.862; 
00256   wc_[2].mean=227.2; wc_[2].sigma=5.070; 
00257   wc_[3].mean=235.7; wc_[3].sigma=6.090; 
00258   wc_[4].mean=243.3; wc_[4].sigma=7.804; 
00259   wc_[5].mean=230.3; wc_[5].sigma=28.91; 
00260 
00261   wc_[6].mean=225.0; wc_[6].sigma=6.000; 
00262   wc_[7].mean=225.0; wc_[7].sigma=6.000; 
00263   wc_[8].mean=225.0; wc_[8].sigma=6.000; 
00264   wc_[9].mean=225.0; wc_[9].sigma=6.000; 
00265 }
00266 */
00267 
00268 void HcalTBTDCUnpacker::reconstructWC(const std::vector<Hit>& hits, HcalTBEventPosition& ep) const {
00269   // process all planes, looping over all hits...
00270   const int MAX_HITS=100;
00271   float hits1[MAX_HITS], hits2[MAX_HITS], hitsA[MAX_HITS];
00272   int n1,n2,nA,chan1,chan2,chanA;
00273   
00274   std::vector<double> x;
00275 
00276   for (int plane=0; plane<PLANECOUNT; plane++) {
00277     n1=0; n2=0; nA=0;
00278 
00279     std::vector<double> plane_hits;
00280     double hit_time;
00281 
00282     chan1=WC_CHANNELIDS[plane*3];
00283     chan2=WC_CHANNELIDS[plane*3+1];
00284     chanA=WC_CHANNELIDS[plane*3+2];
00285 
00286     for (std::vector<Hit>::const_iterator j=hits.begin(); j!=hits.end(); j++) {
00287       if (j->channel==chan1 && n1<MAX_HITS) {
00288         hits1[n1]=j->time-TDC_OFFSET_CONSTANT; n1++;
00289       }
00290       if (j->channel==chan2 && n2<MAX_HITS) {
00291         hits2[n2]=j->time-TDC_OFFSET_CONSTANT; n2++;
00292       }
00293       if (j->channel==chanA && nA<MAX_HITS) {
00294         hitsA[nA]=j->time-TDC_OFFSET_CONSTANT; nA++;
00295       }
00296     }
00297     
00298     // anode-matched hits
00299     for (int ii=0; ii<n1; ii++) {
00300       int jmin=-1, lmin=-1;
00301       float dsumMin=99999;
00302       for (int jj=0; jj<n2; jj++) {
00303         for (int ll=0; ll<nA; ll++) {
00304           float dsum=fabs(wc_[plane].mean - hits1[ii] - hits2[jj] + 2.0*hitsA[ll]);
00305           if(dsum<(N_SIGMA*wc_[plane].sigma) && dsum<dsumMin){
00306             jmin=jj;
00307             lmin=ll;
00308             dsumMin=dsum;
00309            }
00310         }             
00311       }
00312       if (jmin>=0) {
00313         hit_time = wc_[plane].b0 +
00314                              wc_[plane].b1 * (hits1[ii]-hits2[jmin]);
00315         if((plane%2)==0)
00316                 {
00317                 plane_hits.push_back(-hit_time);
00318                 }else{
00319                 plane_hits.push_back(hit_time);
00320                 }
00321         hits1[ii]=-99999;
00322         hits2[jmin]=-99999;
00323         hitsA[lmin]=99999;
00324       }
00325     }
00326 
00327     if (includeUnmatchedHits_||plane>9)   // unmatched hits (all pairs get in here)
00328       for (int ii=0; ii<n1; ii++) {
00329         if (hits1[ii]<-99990) continue;
00330         for (int jj=0; jj<n2; jj++) {
00331           if (hits2[jj]<-99990) continue;
00332           hit_time = wc_[plane].b0 +
00333                              wc_[plane].b1 * (hits1[ii]-hits2[jj]);
00334           if((plane%2)==0)
00335                 {
00336                 plane_hits.push_back(-hit_time);
00337                 }else{
00338                 plane_hits.push_back(hit_time);
00339                 }
00340         }
00341       }
00342 
00343     if ((plane%2)==0) x=plane_hits;
00344     else {
00345       char chamber='A'+plane/2;
00346       ep.setChamberHits(chamber,x,plane_hits);
00347     }
00348   }
00349   
00350 }
00351 
00352 }