Changeset 17862 for trunk/Mars/msignal
- Timestamp:
- 05/15/14 19:09:49 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/Mars/msignal/MTreatSaturation.cc
r17858 r17862 111 111 // Find maximum 112 112 const Float_t *pmax = pbeg; 113 114 113 for (const Float_t *p=pbeg; p<pend; p++) 115 114 if (*p>*pmax) 116 115 pmax = p; 117 116 117 // Is there any chance for saturation? 118 118 if (*pmax<1800) 119 119 continue; 120 120 121 MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray()); 122 s.SetExtractionType(MExtralgoSpline::kAmplitudeRel); 123 124 const double leading = s.SearchYdn(1800., pmax-pbeg); 125 const double falling = s.SearchYup(1800., pmax-pbeg); 126 const double width = falling-leading; 127 const double amplitude = 1800/(0.898417 - 0.0187633*width + 0.000163919*width*width - 6.87369e-7*width*width*width + 1.13264e-9*width*width*width*width); 128 const double time = leading-(4.46944 - 0.277272*width + 0.0035479*width*width); 129 121 // Determine a rough estimate for the average baseline 130 122 double baseline = 0; 131 123 for (const Float_t *p=ptr+10; p<ptr+50; p++) … … 133 125 baseline /= 40; 134 126 135 const UInt_t beg = TMath::CeilNint(leading); 136 const UInt_t end = TMath::FloorNint(falling); 127 // Do a spline interpolation to find the crossing of 1.8V 128 // before and after the maximum 129 MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray()); 130 s.SetExtractionType(MExtralgoSpline::kAmplitudeRel); 137 131 138 ptr += fHiGainFirst+1; 139 for (UInt_t i=0; i<end-beg; i++) 132 const double leading = s.SearchYdn(pmax-pbeg, 1800); 133 const double falling = s.SearchYup(pmax-pbeg, 1800); 134 const double width = falling-leading; 135 136 // If the width above 1.8V is below 3 samples... go ahead as usual. 137 if (width>3) 140 138 { 141 cout << i << " " << ptr[beg-i] << " "; 142 ptr[beg+i] = amplitude*(1-1/(1+exp((i-time)/2.14)))*exp((time-i)/38.8)+baseline; 143 cout << ptr[beg+i] << endl; 139 // Estimate the amplitude and the arrival time from the width 140 const double amplitude = (1800-baseline)/(0.898417 - 0.0187633*width + 0.000163919*width*width - 6.87369e-7*width*width*width + 1.13264e-9*width*width*width*width); 141 const double deltat = -1.41371-0.0525846*width + 93.2763/(width+13.196); 142 const double time0 = leading-deltat-1; 143 144 const Int_t beg = TMath::FloorNint(leading); 145 const Int_t end = TMath::CeilNint(falling); 146 147 ptr += fHiGainFirst+1; 148 for (Int_t i=beg-5; i<end+5; i++) 149 { 150 const double x = i-time0; 151 const double v = amplitude*(1-1/(1+exp(x/2.14)))*exp(-x/38.8)+baseline; 152 if (v>ptr[i]) 153 ptr[i] = v; 154 } 144 155 } 145 146 cout << endl;147 156 } 148 157
Note:
See TracChangeset
for help on using the changeset viewer.