EcalTBHodoscopeReconstructor/src/EcalTBHodoscopeRecInfoAlgo.cc

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