Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:27:15

0001 #include "RecoTBCalo/EcalTBHodoscopeReconstructor/interface/EcalTBHodoscopeRecInfoAlgo.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 
0004 #include <list>
0005 
0006 EcalTBHodoscopeRecInfoAlgo::EcalTBHodoscopeRecInfoAlgo(int fitMethod,
0007                                                        const std::vector<double>& planeShift,
0008                                                        const std::vector<double>& zPosition)
0009     : fitMethod_(fitMethod), planeShift_(planeShift), zPosition_(zPosition) {
0010   myGeometry_ = new EcalTBHodoscopeGeometry();
0011 }
0012 
0013 void EcalTBHodoscopeRecInfoAlgo::clusterPos(
0014     float& x, float& xQuality, const int& ipl, const int& xclus, const int& width) const {
0015   if (width == 1) {
0016     // Single fiber
0017     x = (myGeometry_->getFibreLp(ipl, xclus) + myGeometry_->getFibreRp(ipl, xclus)) / 2.0 - planeShift_[ipl];
0018     xQuality = (myGeometry_->getFibreRp(ipl, xclus) - myGeometry_->getFibreLp(ipl, xclus));
0019   } else if (width == 2) {
0020     // Two half overlapped fibers
0021     x = (myGeometry_->getFibreLp(ipl, xclus + 1) + myGeometry_->getFibreRp(ipl, xclus)) / 2.0 - planeShift_[ipl];
0022     xQuality = (myGeometry_->getFibreRp(ipl, xclus) - myGeometry_->getFibreLp(ipl, xclus + 1));
0023   } else {
0024     // More then two fibers case
0025     x = (myGeometry_->getFibreLp(ipl, xclus) + myGeometry_->getFibreRp(ipl, xclus + width - 1)) / 2.0 -
0026         planeShift_[ipl];
0027     xQuality = (myGeometry_->getFibreRp(ipl, xclus + width - 1) - myGeometry_->getFibreLp(ipl, xclus));
0028   }
0029 }
0030 
0031 void EcalTBHodoscopeRecInfoAlgo::fitHodo(float& x,
0032                                          float& xQuality,
0033                                          const int& ipl,
0034                                          const int& nclus,
0035                                          const std::vector<int>& xclus,
0036                                          const std::vector<int>& wclus) const {
0037   if (nclus == 1) {
0038     // Fill real x as mean position inside the cluster
0039     // Quality - width of cluster
0040     // To calculate sigma one can do sigma=sqrt(Quality**2/12.0)
0041     clusterPos(x, xQuality, ipl, xclus[0], wclus[0]);
0042   } else {
0043     xQuality = -10 - nclus;
0044   }
0045 }
0046 
0047 void EcalTBHodoscopeRecInfoAlgo::fitLine(float& x,
0048                                          float& xSlope,
0049                                          float& xQuality,
0050                                          const int& ipl1,
0051                                          const int& nclus1,
0052                                          const std::vector<int>& xclus1,
0053                                          const std::vector<int>& wclus1,
0054                                          const int& ipl2,
0055                                          const int& nclus2,
0056                                          const std::vector<int>& xclus2,
0057                                          const std::vector<int>& wclus2) const {
0058   if (nclus1 == 0) {  // Fit with one plane
0059     fitHodo(x, xQuality, ipl2, nclus2, xclus2, wclus2);
0060     xSlope = 0.0;  //?? Should we put another number indicating that is not fitted
0061     return;
0062   }
0063   if (nclus2 == 0) {  // Fit with one plane
0064     fitHodo(x, xQuality, ipl1, nclus1, xclus1, wclus1);
0065     xSlope = 0.0;  //?? Should we put another number indicating that is not fitted
0066     return;
0067   }
0068 
0069   // We have clusters in both planes
0070 
0071   float x1, x2, xQ1, xQ2;
0072   float xs, x0, xq;
0073 
0074   std::list<BeamTrack> tracks;
0075 
0076   for (int i1 = 0; i1 < nclus1; i1++) {
0077     for (int i2 = 0; i2 < nclus2; i2++) {
0078       clusterPos(x1, xQ1, ipl1, xclus1[i1], wclus1[i1]);
0079       clusterPos(x2, xQ2, ipl2, xclus2[i2], wclus2[i2]);
0080 
0081       xs = (x2 - x1) / (zPosition_[ipl2] - zPosition_[ipl1]);               // slope
0082       x0 = ((x2 + x1) - xs * (zPosition_[ipl2] + zPosition_[ipl1])) / 2.0;  // x0
0083       xq = (xQ1 + xQ2) / 2.0;                                               // Quality, how i can do better ?
0084       tracks.push_back(BeamTrack(x0, xs, xq));
0085     }
0086   }
0087 
0088   // find track with minimal slope
0089   tracks.sort();
0090 
0091   // Return results
0092   x = tracks.begin()->x;
0093   xSlope = tracks.begin()->xS;
0094   xQuality = tracks.begin()->xQ;
0095 }
0096 
0097 EcalTBHodoscopeRecInfo EcalTBHodoscopeRecInfoAlgo::reconstruct(const EcalTBHodoscopeRawInfo& hodoscopeRawInfo) const {
0098   // Reset Hodo data
0099   float x, y = -100.0;
0100   float xSlope, ySlope = 0.0;
0101   float xQuality, yQuality = -100.0;
0102 
0103   int nclus[4];
0104   std::vector<int> xclus[4];
0105   std::vector<int> wclus[4];
0106 
0107   for (int ipl = 0; ipl < myGeometry_->getNPlanes(); ipl++) {
0108     int nhits = hodoscopeRawInfo[ipl].numberOfFiredHits();
0109     // Finding clusters
0110     nclus[ipl] = 0;
0111     if (nhits > 0) {
0112       int nh = nhits;
0113       int first = 0;
0114       int last = 0;
0115       while (nh > 0) {
0116         while (hodoscopeRawInfo[ipl][first] == 0)
0117           first++;  // start
0118         last = first + 1;
0119         nh--;
0120         do {
0121           while (last < myGeometry_->getNFibres() && hodoscopeRawInfo[ipl][last]) {
0122             last++;
0123             nh--;
0124           }                                                                               //end
0125           if (last + 1 < myGeometry_->getNFibres() && hodoscopeRawInfo[ipl][last + 1]) {  //Skip 1 fibre hole
0126             last += 2;
0127             nh--;
0128             //std::cout << "Skip fibre " << ipl << " " << first << " "<< last << std::endl;
0129           } else {
0130             break;
0131           }
0132         } while (nh > 0 && last < myGeometry_->getNFibres());
0133         wclus[ipl].push_back(last - first);
0134         xclus[ipl].push_back(first);  // Left edge !!!
0135         nclus[ipl]++;
0136 
0137         first = last + 1;
0138       }
0139     }
0140     //    printClusters(ipl);
0141   }
0142 
0143   //! Fit 0
0144   // Straight line fit for one axis
0145   if (fitMethod_ == 0) {
0146     fitLine(x,
0147             xSlope,
0148             xQuality,
0149             0,
0150             nclus[0],
0151             xclus[0],
0152             wclus[0],  // X1
0153             2,
0154             nclus[2],
0155             xclus[2],
0156             wclus[2]);  // X2
0157     fitLine(y,
0158             ySlope,
0159             yQuality,
0160             1,
0161             nclus[1],
0162             xclus[1],
0163             wclus[1],  // Y1
0164             3,
0165             nclus[3],
0166             xclus[3],
0167             wclus[3]);  // Y2
0168   } else if (fitMethod_ == 1) {
0169     //! Fit 1
0170     // x1 and y2 hodoscope
0171     fitHodo(x, xQuality, 0, nclus[0], xclus[0], wclus[0]);  // X1
0172     //   if ( xQuality[1] < 0.0 ) {
0173     //     printFibres(0);
0174     //     printClusters(0);
0175     //   }
0176     fitHodo(y, yQuality, 1, nclus[1], xclus[1], wclus[1]);  // Y1
0177     //   if ( yQuality[1] < 0.0 ) {
0178     //     printFibres(1);
0179     //     printClusters(1);
0180     //   }
0181   } else if (fitMethod_ == 2) {
0182     //! Fit 2
0183     //x2 and y2 hodoscope
0184     fitHodo(x, xQuality, 2, nclus[2], xclus[2], wclus[2]);  // X2
0185     //   if ( xQuality[2] < 0.0 ) {
0186     //     printFibres(2);
0187     //     printClusters(2);
0188     //   }
0189     fitHodo(y, yQuality, 3, nclus[3], xclus[3], wclus[3]);  // Y2
0190     //   if ( yQuality[2] < 0.0 ) {
0191     //     printFibres(3);
0192     //     printClusters(3);
0193     //   }
0194   }
0195 
0196   return EcalTBHodoscopeRecInfo(x, y, xSlope, ySlope, xQuality, yQuality);
0197 }