Changeset 14949 for trunk/Mars/mcore
- Timestamp:
- 02/23/13 12:34:41 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/mcore/DrsCalib.h
r14935 r14949 347 347 // std::cout << " " << it->avg << std::endl; 348 348 349 const Int_t skip = list.size()/10;349 const size_t skip = list.size()/10; 350 350 rc = AverageSteps(list.begin()+skip, list.begin()+list.size()-skip); 351 351 … … 727 727 const float &v1 = val[rel+1];//-avg; 728 728 729 // I s rising edge?729 // If edge is positive ignore all falling edges 730 730 if (edge>0 && v0>0) 731 731 continue; 732 732 733 // I s falling edge?733 // If edge is negative ignore all falling edges 734 734 if (edge<0 && v0<0) 735 735 continue; 736 736 737 // Has sign changed?737 // Check if there is a zero crossing 738 738 if ((v0<0 && v1<0) || (v0>0 && v1>0)) 739 739 continue; 740 740 741 // Calculate the position p of the zero-crossing 742 // within the interval [rel, rel+1] relative to rel 743 // by linear interpolation. 741 744 const double p = v0==v1 ? 0.5 : v0/(v0-v1); 742 745 743 const double l = i+p - (i_prev+p_prev); 744 746 // If this was at least the second zero-crossing detected 745 747 if (i_prev>=0) 746 748 { 749 // Calculate the distance l between the 750 // current and the last zero-crossing 751 const double l = i+p - (i_prev+p_prev); 752 753 // By summation, the average length of each 754 // cell is calculated. For the first and last 755 // fraction of a cell, the fraction is applied 756 // as a weight. 747 757 const double w0 = 1-p_prev; 748 const double w1 = p;749 750 758 fStat[pos+(spos+i_prev)%1024].first += w0*l; 751 759 fStat[pos+(spos+i_prev)%1024].second += w0; … … 757 765 } 758 766 767 const double w1 = p; 759 768 fStat[pos+(spos+i)%1024].first += w1*l; 760 769 fStat[pos+(spos+i)%1024].second += w1; 761 770 } 762 771 772 // Remember this zero-crossing position 763 773 p_prev = p; 764 774 i_prev = i; … … 819 829 const auto end = beg + 1024; 820 830 821 // Calculate mean 831 // First calculate the average length s of a single 832 // zero-crossing interval in the whole range [0;1023] 833 // (which is identical to the/ wavelength of the 834 // calibration signal) 822 835 double s = 0; 823 836 double w = 0; 824 825 837 for (auto it=beg; it!=end; it++) 826 838 { … … 828 840 w += it->second; 829 841 } 830 831 842 s /= w; 832 843 844 // Dividing the average length s of the zero-crossing 845 // interval in the range [0;1023] by the average length 846 // in the interval [0;n] yields the relative size of 847 // the interval in the range [0;n]. 848 // 849 // Example: 850 // Average [0;1023]: 10.00 (global interval size in samples) 851 // Average [0;512]: 8.00 (local interval size in samples) 852 // 853 // Globally, on average one interval is sampled by 10 samples. 854 // In the sub-range [0;512] one interval is sampled on average 855 // by 8 samples. 856 // That means that the interval contains 64 periods, while 857 // in the ideal case (each sample has the same length), it 858 // should contain 51.2 periods. 859 // So, the sampling position 512 corresponds to a time 640, 860 // while in the ideal case with equally spaces samples, 861 // it would correspond to a time 512. 862 // 863 // The offset (defined as 'ideal - real') is then calculated 864 // as 512*(1-10/8) = -128, so that the time is calculated as 865 // 'sampling position minus offset' 866 // 833 867 double sumw = 0; 834 868 double sumv = 0; … … 837 871 // Sums about many values are numerically less stable than 838 872 // just sums over less. So we do the exercise from both sides. 873 // First from the left 839 874 for (auto it=beg; it!=end-512; it++, n++) 840 875 { … … 842 877 const double valw = it->second; 843 878 844 it->first = sumv>0 ? n*(1-s*sumw/sumv) : 0;879 it->first = sumv>0 ? n*(1-s*sumw/sumv) : 0; 845 880 846 881 sumv += valv; … … 852 887 n = 1; 853 888 889 // Second from the right 854 890 for (auto it=end-1; it!=beg-1+512; it--, n++) 855 891 { … … 862 898 it->first = sumv>0 ? n*(s*sumw/sumv-1) : 0; 863 899 } 900 901 // A crosscheck has shown, that the values from the left 902 // and right perfectly agree over the whole range. This means 903 // the a calculation from just one side would be enough, but 904 // doing it from both sides might still make the numerics 905 // a bit more stable. 864 906 } 865 907 }
Note:
See TracChangeset
for help on using the changeset viewer.