00001 // @(#)root/spectrum:$Id: TSpectrum2.cxx 34991 2010-08-25 10:32:26Z brun $ 00002 // Author: Miroslav Morhac 17/01/2006 00003 00004 ///////////////////////////////////////////////////////////////////////////// 00005 // THIS CLASS CONTAINS ADVANCED SPECTRA PROCESSING FUNCTIONS. // 00006 // // 00007 // ONE-DIMENSIONAL BACKGROUND ESTIMATION FUNCTIONS // 00008 // TWO-DIMENSIONAL BACKGROUND ESTIMATION FUNCTIONS // 00009 // ONE-DIMENSIONAL SMOOTHING FUNCTIONS // 00010 // TWO-DIMENSIONAL SMOOTHING FUNCTIONS // 00011 // ONE-DIMENSIONAL DECONVOLUTION FUNCTIONS // 00012 // TWO-DIMENSIONAL DECONVOLUTION FUNCTIONS // 00013 // ONE-DIMENSIONAL PEAK SEARCH FUNCTIONS // 00014 // TWO-DIMENSIONAL PEAK SEARCH FUNCTIONS // 00015 // // 00016 // These functions were written by: // 00017 // Miroslav Morhac // 00018 // Institute of Physics // 00019 // Slovak Academy of Sciences // 00020 // Dubravska cesta 9, 842 28 BRATISLAVA // 00021 // SLOVAKIA // 00022 // // 00023 // email:fyzimiro@savba.sk, fax:+421 7 54772479 // 00024 // // 00025 // The original code in C has been repackaged as a C++ class by R.Brun // 00026 // // 00027 // The algorithms in this class have been published in the following // 00028 // references: // 00029 // [1] M.Morhac et al.: Background elimination methods for // 00030 // multidimensional coincidence gamma-ray spectra. Nuclear // 00031 // Instruments and Methods in Physics Research A 401 (1997) 113- // 00032 // 132. // 00033 // // 00034 // [2] M.Morhac et al.: Efficient one- and two-dimensional Gold // 00035 // deconvolution and its application to gamma-ray spectra // 00036 // decomposition. Nuclear Instruments and Methods in Physics // 00037 // Research A 401 (1997) 385-408. // 00038 // // 00039 // [3] M.Morhac et al.: Identification of peaks in multidimensional // 00040 // coincidence gamma-ray spectra. Nuclear Instruments and Methods in // 00041 // Research Physics A 443(2000), 108-125. // 00042 // // 00043 // These NIM papers are also available as Postscript files from: // 00044 // 00045 /* 00046 ftp://root.cern.ch/root/SpectrumDec.ps.gz 00047 ftp://root.cern.ch/root/SpectrumSrc.ps.gz 00048 ftp://root.cern.ch/root/SpectrumBck.ps.gz 00049 */ 00050 // 00051 ///////////////////////////////////////////////////////////////////////////// 00052 // 00053 /////////////////////NEW FUNCTIONS January 2006 00054 //Begin_Html <!-- 00055 /* --> 00056 <div class=Section1> 00057 00058 <p class=MsoNormal style='text-align:justify'><i><span lang=EN-US 00059 style='font-size:16.0pt'>All figures in this page were prepared using DaqProVis 00060 system, Data Acquisition, Processing and Visualization system, which is being 00061 developed at the Institute of Physics, Slovak Academy of Sciences, Bratislava, 00062 Slovakia: </span></i></p> 00063 00064 <p class=MsoNormal style='text-align:justify'><i><span lang=EN-US 00065 style='font-size:16.0pt'><a href="http://www.fu.sav.sk/nph/projects/DaqProVis/">http://www.fu.sav.sk/nph/projects/DaqProVis/</a> 00066 under construction</span></i></p> 00067 00068 <p class=MsoNormal style='text-align:justify'><i><span lang=EN-US 00069 style='font-size:16.0pt'><a href="http://www.fu.sav.sk/nph/projects/ProcFunc/">http://www.fu.sav.sk/nph/projects/ProcFunc/</a> 00070 .</span></i></p> 00071 00072 </div> 00073 00074 <!-- */ 00075 // --> End_Html 00076 00077 //______________________________________________________________________________ 00078 00079 00080 00081 #include "TSpectrum2.h" 00082 #include "TPolyMarker.h" 00083 #include "TList.h" 00084 #include "TH1.h" 00085 #include "TMath.h" 00086 #define PEAK_WINDOW 1024 00087 00088 Int_t TSpectrum2::fgIterations = 3; 00089 Int_t TSpectrum2::fgAverageWindow = 3; 00090 00091 ClassImp(TSpectrum2) 00092 00093 //______________________________________________________________________________ 00094 TSpectrum2::TSpectrum2() :TNamed("Spectrum", "Miroslav Morhac peak finder") 00095 { 00096 // Constructor. 00097 00098 Int_t n = 100; 00099 fMaxPeaks = n; 00100 fPosition = new Float_t[n]; 00101 fPositionX = new Float_t[n]; 00102 fPositionY = new Float_t[n]; 00103 fResolution = 1; 00104 fHistogram = 0; 00105 fNPeaks = 0; 00106 } 00107 00108 00109 //______________________________________________________________________________ 00110 TSpectrum2::TSpectrum2(Int_t maxpositions, Float_t resolution) :TNamed("Spectrum", "Miroslav Morhac peak finder") 00111 { 00112 // maxpositions: maximum number of peaks 00113 // resolution: determines resolution of the neighboring peaks 00114 // default value is 1 correspond to 3 sigma distance 00115 // between peaks. Higher values allow higher resolution 00116 // (smaller distance between peaks. 00117 // May be set later through SetResolution. 00118 00119 Int_t n = maxpositions; 00120 fMaxPeaks = n; 00121 fPosition = new Float_t[n]; 00122 fPositionX = new Float_t[n]; 00123 fPositionY = new Float_t[n]; 00124 fHistogram = 0; 00125 fNPeaks = 0; 00126 SetResolution(resolution); 00127 } 00128 00129 00130 //______________________________________________________________________________ 00131 TSpectrum2::~TSpectrum2() 00132 { 00133 // Destructor. 00134 00135 delete [] fPosition; 00136 delete [] fPositionX; 00137 delete [] fPositionY; 00138 delete fHistogram; 00139 } 00140 00141 00142 //______________________________________________________________________________ 00143 void TSpectrum2::SetAverageWindow(Int_t w) 00144 { 00145 // static function: Set average window of searched peaks 00146 // see TSpectrum2::SearchHighRes 00147 00148 fgAverageWindow = w; 00149 } 00150 00151 //______________________________________________________________________________ 00152 void TSpectrum2::SetDeconIterations(Int_t n) 00153 { 00154 // static function: Set max number of decon iterations in deconvolution operation 00155 // see TSpectrum2::SearchHighRes 00156 00157 fgIterations = n; 00158 } 00159 00160 00161 //______________________________________________________________________________ 00162 TH1 *TSpectrum2::Background(const TH1 * h, int number_of_iterations, 00163 Option_t * option) 00164 { 00165 ///////////////////////////////////////////////////////////////////////////// 00166 // TWO-DIMENSIONAL BACKGROUND ESTIMATION FUNCTION // 00167 // This function calculates the background spectrum in the input histogram h. 00168 // The background is returned as a histogram. 00169 // 00170 // Function parameters: 00171 // -h: input 2-d histogram 00172 // -numberIterations, (default value = 20) 00173 // Increasing numberIterations make the result smoother and lower. 00174 // -option: may contain one of the following options 00175 // - to set the direction parameter 00176 // "BackIncreasingWindow". By default the direction is BackDecreasingWindow 00177 // - filterOrder-order of clipping filter, (default "BackOrder2" 00178 // -possible values= "BackOrder4" 00179 // "BackOrder6" 00180 // "BackOrder8" 00181 // - "nosmoothing"- if selected, the background is not smoothed 00182 // By default the background is smoothed. 00183 // - smoothWindow-width of smoothing window, (default is "BackSmoothing3") 00184 // -possible values= "BackSmoothing5" 00185 // "BackSmoothing7" 00186 // "BackSmoothing9" 00187 // "BackSmoothing11" 00188 // "BackSmoothing13" 00189 // "BackSmoothing15" 00190 // - "Compton" if selected the estimation of Compton edge 00191 // will be included. 00192 // - "same" : if this option is specified, the resulting background 00193 // histogram is superimposed on the picture in the current pad. 00194 // 00195 // NOTE that the background is only evaluated in the current range of h. 00196 // ie, if h has a bin range (set via h->GetXaxis()->SetRange(binmin,binmax), 00197 // the returned histogram will be created with the same number of bins 00198 // as the input histogram h, but only bins from binmin to binmax will be filled 00199 // with the estimated background. 00200 // // 00201 ///////////////////////////////////////////////////////////////////////////// 00202 Error("Background","function not yet implemented: h=%s, iter=%d, option=%sn" 00203 , h->GetName(), number_of_iterations, option); 00204 return 0; 00205 } 00206 00207 //______________________________________________________________________________ 00208 void TSpectrum2::Print(Option_t *) const 00209 { 00210 // Print the array of positions 00211 00212 printf("\nNumber of positions = %d\n",fNPeaks); 00213 for (Int_t i=0;i<fNPeaks;i++) { 00214 printf(" x[%d] = %g, y[%d] = %g\n",i,fPositionX[i],i,fPositionY[i]); 00215 } 00216 } 00217 00218 00219 00220 //______________________________________________________________________________ 00221 Int_t TSpectrum2::Search(const TH1 * hin, Double_t sigma, 00222 Option_t * option, Double_t threshold) 00223 { 00224 ///////////////////////////////////////////////////////////////////////////// 00225 // TWO-DIMENSIONAL PEAK SEARCH FUNCTION // 00226 // This function searches for peaks in source spectrum in hin // 00227 // The number of found peaks and their positions are written into // 00228 // the members fNpeaks and fPositionX. // 00229 // The search is performed in the current histogram range. // 00230 // // 00231 // Function parameters: // 00232 // hin: pointer to the histogram of source spectrum // 00233 // sigma: sigma of searched peaks, for details we refer to manual // 00234 // threshold: (default=0.05) peaks with amplitude less than // 00235 // threshold*highest_peak are discarded. 0<threshold<1 // 00236 // // 00237 // By default, the background is removed before deconvolution. // 00238 // Specify the option "nobackground" to not remove the background. // // 00239 // // 00240 // By default the "Markov" chain algorithm is used. // 00241 // Specify the option "noMarkov" to disable this algorithm // 00242 // Note that by default the source spectrum is replaced by a new spectrum// // 00243 // // 00244 // By default a polymarker object is created and added to the list of // 00245 // functions of the histogram. The histogram is drawn with the specified // 00246 // option and the polymarker object drawn on top of the histogram. // 00247 // The polymarker coordinates correspond to the npeaks peaks found in // 00248 // the histogram. // 00249 // A pointer to the polymarker object can be retrieved later via: // 00250 // TList *functions = hin->GetListOfFunctions(); // 00251 // TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker") // 00252 // Specify the option "goff" to disable the storage and drawing of the // 00253 // polymarker. // 00254 // // 00255 ///////////////////////////////////////////////////////////////////////////// 00256 00257 if (hin == 0) 00258 return 0; 00259 Int_t dimension = hin->GetDimension(); 00260 if (dimension != 2) { 00261 Error("Search", "Must be a 2-d histogram"); 00262 return 0; 00263 } 00264 00265 TString opt = option; 00266 opt.ToLower(); 00267 Bool_t background = kTRUE; 00268 if (opt.Contains("nobackground")) { 00269 background = kFALSE; 00270 opt.ReplaceAll("nobackground",""); 00271 } 00272 Bool_t markov = kTRUE; 00273 if (opt.Contains("nomarkov")) { 00274 markov = kFALSE; 00275 opt.ReplaceAll("nomarkov",""); 00276 } 00277 00278 Int_t sizex = hin->GetXaxis()->GetNbins(); 00279 Int_t sizey = hin->GetYaxis()->GetNbins(); 00280 Int_t i, j, binx,biny, npeaks; 00281 Float_t ** source = new float *[sizex]; 00282 Float_t ** dest = new float *[sizex]; 00283 for (i = 0; i < sizex; i++) { 00284 source[i] = new float[sizey]; 00285 dest[i] = new float[sizey]; 00286 for (j = 0; j < sizey; j++) { 00287 source[i][j] = (Float_t) hin->GetBinContent(i + 1, j + 1); 00288 } 00289 } 00290 //npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, kTRUE, 3, kTRUE, 10); 00291 //the smoothing option is used for 1-d but not for 2-d histograms 00292 npeaks = SearchHighRes(source, dest, sizex, sizey, sigma, 100*threshold, background, fgIterations, markov, fgAverageWindow); 00293 00294 //The logic in the loop should be improved to use the fact 00295 //that fPositionX,Y give a precise position inside a bin. 00296 //The current algorithm takes the center of the bin only. 00297 for (i = 0; i < npeaks; i++) { 00298 binx = 1 + Int_t(fPositionX[i] + 0.5); 00299 biny = 1 + Int_t(fPositionY[i] + 0.5); 00300 fPositionX[i] = hin->GetXaxis()->GetBinCenter(binx); 00301 fPositionY[i] = hin->GetYaxis()->GetBinCenter(biny); 00302 } 00303 for (i = 0; i < sizex; i++) { 00304 delete [] source[i]; 00305 delete [] dest[i]; 00306 } 00307 delete [] source; 00308 delete [] dest; 00309 00310 if (opt.Contains("goff")) 00311 return npeaks; 00312 if (!npeaks) return 0; 00313 TPolyMarker * pm = (TPolyMarker*)hin->GetListOfFunctions()->FindObject("TPolyMarker"); 00314 if (pm) { 00315 hin->GetListOfFunctions()->Remove(pm); 00316 delete pm; 00317 } 00318 pm = new TPolyMarker(npeaks, fPositionX, fPositionY); 00319 hin->GetListOfFunctions()->Add(pm); 00320 pm->SetMarkerStyle(23); 00321 pm->SetMarkerColor(kRed); 00322 pm->SetMarkerSize(1.3); 00323 ((TH1*)hin)->Draw(option); 00324 return npeaks; 00325 } 00326 00327 00328 //______________________________________________________________________________ 00329 void TSpectrum2::SetResolution(Float_t resolution) 00330 { 00331 // resolution: determines resolution of the neighboring peaks 00332 // default value is 1 correspond to 3 sigma distance 00333 // between peaks. Higher values allow higher resolution 00334 // (smaller distance between peaks. 00335 // May be set later through SetResolution. 00336 if (resolution > 1) 00337 fResolution = resolution; 00338 else 00339 fResolution = 1; 00340 } 00341 00342 00343 //_____________________________________________________________________________ 00344 //_____________________________________________________________________________ 00345 00346 /////////////////////NEW FUNCTIONS JANUARY 2006 00347 //______________________________________________________________________________ 00348 const char *TSpectrum2::Background(float **spectrum, 00349 Int_t ssizex, Int_t ssizey, 00350 Int_t numberIterationsX, 00351 Int_t numberIterationsY, 00352 Int_t direction, 00353 Int_t filterType) 00354 { 00355 ///////////////////////////////////////////////////////////////////////////// 00356 // TWO-DIMENSIONAL BACKGROUND ESTIMATION FUNCTION - RECTANGULAR RIDGES // 00357 // This function calculates background spectrum from source spectrum. // 00358 // The result is placed to the array pointed by spectrum pointer. // 00359 // // 00360 // Function parameters: // 00361 // spectrum-pointer to the array of source spectrum // 00362 // ssizex-x length of spectrum // 00363 // ssizey-y length of spectrum // 00364 // numberIterationsX-maximal x width of clipping window // 00365 // numberIterationsY-maximal y width of clipping window // 00366 // for details we refer to manual // 00367 // direction- direction of change of clipping window // 00368 // - possible values=kBackIncreasingWindow // 00369 // kBackDecreasingWindow // 00370 // filterType-determines the algorithm of the filtering // 00371 // -possible values=kBackSuccessiveFiltering // 00372 // kBackOneStepFiltering // 00373 // // 00374 // // 00375 ///////////////////////////////////////////////////////////////////////////// 00376 // 00377 //Begin_Html <!-- 00378 /* --> 00379 <div class=Section1> 00380 00381 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:20.0pt'>Background 00382 estimation</span></b></p> 00383 00384 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p> 00385 00386 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>Goal: 00387 Separation of useful information (peaks) from useless information (background)</span></i><span 00388 style='font-size:18.0pt'> </span></p> 00389 00390 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 00391 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 00392 </span></span><span style='font-size:16.0pt'>method is based on Sensitive 00393 Nonlinear Iterative Peak (SNIP) clipping algorithm [1]</span></p> 00394 00395 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 00396 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 00397 </span></span><span style='font-size:16.0pt'>there exist two algorithms for the 00398 estimation of new value in the channel “<sub><img width=28 height=24 00399 src="gif/TSpectrum2_Background1.gif"></sub>”</span></p> 00400 00401 <p class=MsoNormal style='margin-left:18.0pt;text-align:justify'><span 00402 style='font-size:16.0pt'> </span></p> 00403 00404 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'>Algorithm 00405 based on Successive Comparisons</span></i></p> 00406 00407 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>It 00408 is an extension of one-dimensional SNIP algorithm to another dimension. For 00409 details we refer to [2].</span></p> 00410 00411 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 00412 00413 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:16.0pt'>Algorithm 00414 based on One Step Filtering</span></i></p> 00415 00416 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>New 00417 value in the estimated channel is calculated as</span></p> 00418 00419 <p class=MsoNormal style='text-align:justify'><sub><img width=133 height=39 00420 src="gif/TSpectrum2_Background2.gif"></sub></p> 00421 00422 <p class=MsoNormal style='text-align:justify'> </p> 00423 00424 <p class=MsoNormal style='text-align:justify'> </p> 00425 00426 <p class=MsoNormal style='text-align:justify'><sub><img width=600 height=128 00427 src="gif/TSpectrum2_Background3.gif"></sub></p> 00428 00429 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 00430 00431 <p class=MsoNormal><sub><img width=190 height=38 00432 src="gif/TSpectrum2_Background4.gif"></sub>.</p> 00433 00434 <p class=MsoNormal> </p> 00435 00436 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>where 00437 p = 1, 2, …, number_of_iterations. </span></p> 00438 00439 <p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p> 00440 00441 <p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p> 00442 00443 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'>const 00444 <a href="http://root.cern.ch/root/html/ListOfTypes.html#char" target="_parent">char</a>* 00445 </span></b><span style='font-size:18.0pt'><a 00446 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrum2::Background</b></a><b> 00447 (<a href="http://root.cern.ch/root/html/ListOfTypes.html#float" target="_parent">float</a> 00448 **spectrum, <a href="http://root.cern.ch/root/html/ListOfTypes.html#int" 00449 target="_parent">int</a> ssizex, <a 00450 href="http://root.cern.ch/root/html/ListOfTypes.html#int" target="_parent">int</a> 00451 ssizey, <a href="http://root.cern.ch/root/html/ListOfTypes.html#int" 00452 target="_parent">int</a> numberIterationsX, <a 00453 href="http://root.cern.ch/root/html/ListOfTypes.html#int" target="_parent">int</a> 00454 numberIterationsY, <a href="http://root.cern.ch/root/html/ListOfTypes.html#int" 00455 target="_parent">int</a> direction, <a 00456 href="http://root.cern.ch/root/html/ListOfTypes.html#int" target="_parent">int</a> 00457 filterType) </b></span></p> 00458 00459 <p class=MsoNormal> </p> 00460 00461 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This 00462 function calculates background spectrum from the source spectrum. The result 00463 is placed in the matrix pointed by spectrum pointer. One can also switch the 00464 direction of the change of the clipping window and to select one of the two 00465 above given algorithms.</span><span style='font-size:18.0pt'> </span><span 00466 style='font-size:16.0pt'>On successful completion it returns 0. On error it 00467 returns pointer to the string describing error.</span></p> 00468 00469 <p class=MsoNormal> </p> 00470 00471 <p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p> 00472 00473 <p class=MsoNormal> <b><span style='font-size:14.0pt'>spectrum</span></b>-pointer 00474 to the matrix of source spectrum </p> 00475 00476 <p class=MsoNormal> <b><span style='font-size:14.0pt'>ssizex, ssizey</span></b>-lengths 00477 of the spectrum matrix </p> 00478 00479 <p class=MsoNormal style='text-align:justify'> <b><span 00480 style='font-size:14.0pt'>numberIterationsX, numberIterationsY</span></b>maximal 00481 widths of clipping</p> 00482 00483 <p class=MsoNormal style='text-align:justify'> window, 00484 </p> 00485 00486 <p class=MsoNormal> <b><span style='font-size:14.0pt'>direction</span></b>- 00487 direction of change of clipping window </p> 00488 00489 <p class=MsoNormal> - possible 00490 values=kBackIncreasingWindow </p> 00491 00492 <p class=MsoNormal> 00493 kBackDecreasingWindow </p> 00494 00495 <p class=MsoNormal> <b><span style='font-size:14.0pt'>filterType</span></b>-type 00496 of the clipping algorithm, </p> 00497 00498 <p class=MsoNormal> -possible values=kBack SuccessiveFiltering</p> 00499 00500 <p class=MsoNormal> 00501 kBackOneStepFiltering </p> 00502 00503 <p class=MsoNormal> </p> 00504 00505 <p class=MsoNormal><b><i><span style='font-size:18.0pt'>References:</span></i></b></p> 00506 00507 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[1] 00508 C. G Ryan et al.: SNIP, a statistics-sensitive background treatment for the 00509 quantitative analysis of PIXE spectra in geoscience applications. NIM, B34 00510 (1988), 396-402.</span></p> 00511 00512 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[2] 00513 </span><span lang=SK style='font-size:16.0pt'> M. Morháč, J. Kliman, V. 00514 Matoušek, M. Veselský, I. Turzo</span><span style='font-size:16.0pt'>.: 00515 Background elimination methods for multidimensional gamma-ray spectra. NIM, 00516 A401 (1997) 113-132.</span></p> 00517 00518 </div> 00519 00520 <!-- */ 00521 // --> End_Html 00522 //Begin_Html <!-- 00523 /* --> 00524 <div class=Section1> 00525 00526 <p class=MsoNormal><i><span style='font-size:18.0pt'>Example 1– script Back_gamma64.c 00527 :</span></i></p> 00528 00529 <p class=MsoNormal><img width=602 height=418 00530 src="gif/TSpectrum2_Background1.jpg"></p> 00531 00532 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00533 1 Original two-dimensional gamma-gamma-ray spectrum</span></b></p> 00534 00535 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 00536 width=602 height=418 src="gif/TSpectrum2_Background2.jpg"></span></b></p> 00537 00538 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00539 2 Background estimated from data from Fig. 1 using decreasing clipping window with 00540 widths 4, 4 and algorithm based on successive comparisons. The estimate 00541 includes not only continuously changing background but also one-dimensional 00542 ridges.</span></b></p> 00543 00544 <p class=MsoNormal><b><span style='font-size:14.0pt;color:green'><img 00545 width=602 height=418 src="gif/TSpectrum2_Background3.jpg"></span></b></p> 00546 00547 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00548 3 Resulting peaks after subtraction of the estimated background (Fig. 2) from original 00549 two-dimensional gamma-gamma-ray spectrum (Fig. 1).</span></b></p> 00550 00551 <p class=MsoNormal><b><span style='font-size:14.0pt;color:green'> </span></b></p> 00552 00553 <p class=MsoNormal><b><span style='font-size:14.0pt;color:green'> </span></b></p> 00554 00555 <p class=MsoNormal><b><span style='font-size:16.0pt;color:green'>Script:</span></b></p> 00556 00557 <p class=MsoNormal>// Example to illustrate the background estimator (class 00558 TSpectrum).</p> 00559 00560 <p class=MsoNormal>// To execute this example, do</p> 00561 00562 <p class=MsoNormal>// root > .x Back_gamma64.C</p> 00563 00564 <p class=MsoNormal> </p> 00565 00566 <p class=MsoNormal>#include <TSpectrum> </p> 00567 00568 <p class=MsoNormal> </p> 00569 00570 <p class=MsoNormal>void Back_gamma64() {</p> 00571 00572 <p class=MsoNormal> Int_t i, j;</p> 00573 00574 <p class=MsoNormal> Double_t nbinsx = 64;</p> 00575 00576 <p class=MsoNormal> Double_t nbinsy = 64; </p> 00577 00578 <p class=MsoNormal> Double_t xmin = 0;</p> 00579 00580 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 00581 00582 <p class=MsoNormal> Double_t ymin = 0;</p> 00583 00584 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 00585 00586 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 00587 00588 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 00589 00590 <p class=MsoNormal> source[i]=new float[nbinsy]; </p> 00591 00592 <p class=MsoNormal> TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 00593 00594 <p class=MsoNormal> TFile *f = new 00595 TFile("spectra2\\TSpectrum2.root");</p> 00596 00597 <p class=MsoNormal> back=(TH2F*) f->Get("back1;1");</p> 00598 00599 <p class=MsoNormal> TCanvas *Background = new 00600 TCanvas("Background","Estimation of background with increasing 00601 window",10,10,1000,700);</p> 00602 00603 <p class=MsoNormal> TSpectrum *s = new TSpectrum();</p> 00604 00605 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 00606 00607 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 00608 00609 <p class=MsoNormal> source[i][j] = back->GetBinContent(i + 00610 1,j + 1); </p> 00611 00612 <p class=MsoNormal> }</p> 00613 00614 <p class=MsoNormal> } </p> 00615 00616 <p class=MsoNormal> s->Background(source,nbinsx,nbinsy,4,4,kBackDecreasingWindow,kBackSuccessiveFiltering);</p> 00617 00618 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 00619 00620 <p class=MsoNormal> for (j = 0; j < nbinsy; j++)</p> 00621 00622 <p class=MsoNormal> back->SetBinContent(i + 1,j + 1, source[i][j]); </p> 00623 00624 <p class=MsoNormal> }</p> 00625 00626 <p class=MsoNormal> back->Draw("SURF"); </p> 00627 00628 <p class=MsoNormal> }</p> 00629 00630 <p class=MsoNormal> </p> 00631 00632 <p class=MsoNormal><i><span style='font-size:18.0pt'>Example 2– script Back_gamma256.c 00633 :</span></i></p> 00634 00635 <p class=MsoNormal><img width=602 height=418 00636 src="gif/TSpectrum2_Background4.jpg"></p> 00637 00638 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00639 4 Original two-dimensional gamma-gamma-ray spectrum 256x256 channels</span></b></p> 00640 00641 <p class=MsoNormal><img width=602 height=418 00642 src="gif/TSpectrum2_Background5.jpg"></p> 00643 00644 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00645 5 Peaks after subtraction of the estimated background (increasing clipping 00646 window with widths 8, 8 and algorithm based on successive comparisons) from original 00647 two-dimensional gamma-gamma-ray spectrum (Fig. 4).</span></b></p> 00648 00649 <p class=MsoNormal><b><span style='font-size:16.0pt;color:green'> </span></b></p> 00650 00651 <p class=MsoNormal><b><span style='font-size:16.0pt;color:green'>Script:</span></b></p> 00652 00653 <p class=MsoNormal>// Example to illustrate the background estimator (class 00654 TSpectrum).</p> 00655 00656 <p class=MsoNormal>// To execute this example, do</p> 00657 00658 <p class=MsoNormal>// root > .x Back_gamma256.C</p> 00659 00660 <p class=MsoNormal> </p> 00661 00662 <p class=MsoNormal>#include <TSpectrum> </p> 00663 00664 <p class=MsoNormal> </p> 00665 00666 <p class=MsoNormal>void Back_gamma256() {</p> 00667 00668 <p class=MsoNormal> Int_t i, j;</p> 00669 00670 <p class=MsoNormal> Double_t nbinsx = 64;</p> 00671 00672 <p class=MsoNormal> Double_t nbinsy = 64; </p> 00673 00674 <p class=MsoNormal> Double_t xmin = 0;</p> 00675 00676 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 00677 00678 <p class=MsoNormal> Double_t ymin = 0;</p> 00679 00680 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 00681 00682 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 00683 00684 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 00685 00686 <p class=MsoNormal> source[i]=new float[nbinsy]; </p> 00687 00688 <p class=MsoNormal> TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 00689 00690 <p class=MsoNormal> TFile *f = new 00691 TFile("spectra2\\TSpectrum2.root");</p> 00692 00693 <p class=MsoNormal> back=(TH2F*) f->Get("back2;1");</p> 00694 00695 <p class=MsoNormal> TCanvas *Background = new 00696 TCanvas("Background","Estimation of background with increasing 00697 window",10,10,1000,700);</p> 00698 00699 <p class=MsoNormal> TSpectrum *s = new TSpectrum();</p> 00700 00701 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 00702 00703 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 00704 00705 <p class=MsoNormal> source[i][j] = back->GetBinContent(i + 00706 1,j + 1); </p> 00707 00708 <p class=MsoNormal> }</p> 00709 00710 <p class=MsoNormal> } </p> 00711 00712 <p class=MsoNormal> s->Background(source,nbinsx,nbinsy,8,8,kBackIncreasingWindow,kBackSuccessiveFiltering);</p> 00713 00714 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 00715 00716 <p class=MsoNormal> for (j = 0; j < nbinsy; j++)</p> 00717 00718 <p class=MsoNormal> back->SetBinContent(i + 1,j + 1, source[i][j]); </p> 00719 00720 <p class=MsoNormal> }</p> 00721 00722 <p class=MsoNormal> back->Draw("SURF"); </p> 00723 00724 <p class=MsoNormal> }</p> 00725 00726 <p class=MsoNormal><i><span style='font-size:18.0pt'>Example 3– script Back_synt256.c 00727 :</span></i></p> 00728 00729 <p class=MsoNormal><img width=602 height=418 00730 src="gif/TSpectrum2_Background6.jpg"></p> 00731 00732 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00733 6 Original two-dimensional synthetic spectrum 256x256 channels</span></b></p> 00734 00735 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 00736 width=602 height=418 src="gif/TSpectrum2_Background7.jpg"></span></b></p> 00737 00738 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00739 7 Peaks after subtraction of the estimated background (increasing clipping 00740 window with widths 8, 8 and algorithm based on successive comparisons) from original 00741 two-dimensional gamma-gamma-ray spectrum (Fig. 6). One can observe artifacts 00742 (ridges) around the peaks due to the employed algorithm.</span></b></p> 00743 00744 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 00745 width=602 height=418 src="gif/TSpectrum2_Background8.jpg"></span></b></p> 00746 00747 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 00748 8 Peaks after subtraction of the estimated background (increasing clipping 00749 window with widths 8, 8 and algorithm based on one step filtering) from original 00750 two-dimensional gamma-gamma-ray spectrum (Fig. 6). The artifacts from the 00751 above given Fig. 7 disappeared.</span></b></p> 00752 00753 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'> </span></b></p> 00754 00755 <p class=MsoNormal><b><span style='font-size:16.0pt;color:green'>Script:</span></b></p> 00756 00757 <p class=MsoNormal>// Example to illustrate the background estimator (class 00758 TSpectrum).</p> 00759 00760 <p class=MsoNormal>// To execute this example, do</p> 00761 00762 <p class=MsoNormal>// root > .x Back_synt256.C</p> 00763 00764 <p class=MsoNormal> </p> 00765 00766 <p class=MsoNormal>#include <TSpectrum> </p> 00767 00768 <p class=MsoNormal> </p> 00769 00770 <p class=MsoNormal>void Back_synt256() {</p> 00771 00772 <p class=MsoNormal> Int_t i, j;</p> 00773 00774 <p class=MsoNormal> Double_t nbinsx = 64;</p> 00775 00776 <p class=MsoNormal> Double_t nbinsy = 64; </p> 00777 00778 <p class=MsoNormal> Double_t xmin = 0;</p> 00779 00780 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 00781 00782 <p class=MsoNormal> Double_t ymin = 0;</p> 00783 00784 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 00785 00786 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 00787 00788 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 00789 00790 <p class=MsoNormal> source[i]=new float[nbinsy]; </p> 00791 00792 <p class=MsoNormal> TH2F *back = new TH2F("back","Background estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 00793 00794 <p class=MsoNormal> TFile *f = new 00795 TFile("spectra2\\TSpectrum2.root");</p> 00796 00797 <p class=MsoNormal> back=(TH2F*) f->Get("back3;1");</p> 00798 00799 <p class=MsoNormal> TCanvas *Background = new 00800 TCanvas("Background","Estimation of background with increasing 00801 window",10,10,1000,700);</p> 00802 00803 <p class=MsoNormal> TSpectrum *s = new TSpectrum();</p> 00804 00805 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 00806 00807 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 00808 00809 <p class=MsoNormal> source[i][j] = back->GetBinContent(i + 00810 1,j + 1); </p> 00811 00812 <p class=MsoNormal> }</p> 00813 00814 <p class=MsoNormal> } </p> 00815 00816 <p class=MsoNormal> s->Background(source,nbinsx,nbinsy,8,8,kBackIncreasingWindow,kBackSuccessiveFiltering);//kBackOneStepFiltering</p> 00817 00818 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 00819 00820 <p class=MsoNormal> for (j = 0; j < nbinsy; j++)</p> 00821 00822 <p class=MsoNormal> back->SetBinContent(i + 1,j + 1, source[i][j]); </p> 00823 00824 <p class=MsoNormal> }</p> 00825 00826 <p class=MsoNormal> back->Draw("SURF"); </p> 00827 00828 <p class=MsoNormal> }</p> 00829 00830 </div> 00831 00832 <!-- */ 00833 // --> End_Html 00834 00835 int i, x, y, sampling, r1, r2; 00836 float a, b, p1, p2, p3, p4, s1, s2, s3, s4; 00837 if (ssizex <= 0 || ssizey <= 0) 00838 return "Wrong parameters"; 00839 if (numberIterationsX < 1 || numberIterationsY < 1) 00840 return "Width of Clipping Window Must Be Positive"; 00841 if (ssizex < 2 * numberIterationsX + 1 00842 || ssizey < 2 * numberIterationsY + 1) 00843 return ("Too Large Clipping Window"); 00844 float **working_space = new float *[ssizex]; 00845 for (i = 0; i < ssizex; i++) 00846 working_space[i] = new float[ssizey]; 00847 sampling = 00848 (int) TMath::Max(numberIterationsX, numberIterationsY); 00849 if (direction == kBackIncreasingWindow) { 00850 if (filterType == kBackSuccessiveFiltering) { 00851 for (i = 1; i <= sampling; i++) { 00852 r1 = (int) TMath::Min(i, numberIterationsX), r2 = 00853 (int) TMath::Min(i, numberIterationsY); 00854 for (y = r2; y < ssizey - r2; y++) { 00855 for (x = r1; x < ssizex - r1; x++) { 00856 a = spectrum[x][y]; 00857 p1 = spectrum[x - r1][y - r2]; 00858 p2 = spectrum[x - r1][y + r2]; 00859 p3 = spectrum[x + r1][y - r2]; 00860 p4 = spectrum[x + r1][y + r2]; 00861 s1 = spectrum[x][y - r2]; 00862 s2 = spectrum[x - r1][y]; 00863 s3 = spectrum[x + r1][y]; 00864 s4 = spectrum[x][y + r2]; 00865 b = (p1 + p2) / 2.0; 00866 if (b > s2) 00867 s2 = b; 00868 b = (p1 + p3) / 2.0; 00869 if (b > s1) 00870 s1 = b; 00871 b = (p2 + p4) / 2.0; 00872 if (b > s4) 00873 s4 = b; 00874 b = (p3 + p4) / 2.0; 00875 if (b > s3) 00876 s3 = b; 00877 s1 = s1 - (p1 + p3) / 2.0; 00878 s2 = s2 - (p1 + p2) / 2.0; 00879 s3 = s3 - (p3 + p4) / 2.0; 00880 s4 = s4 - (p2 + p4) / 2.0; 00881 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + 00882 p3 + 00883 p4) / 4.0; 00884 if (b < a && b > 0) 00885 a = b; 00886 working_space[x][y] = a; 00887 } 00888 } 00889 for (y = r2; y < ssizey - r2; y++) { 00890 for (x = r1; x < ssizex - r1; x++) { 00891 spectrum[x][y] = working_space[x][y]; 00892 } 00893 } 00894 } 00895 } 00896 00897 else if (filterType == kBackOneStepFiltering) { 00898 for (i = 1; i <= sampling; i++) { 00899 r1 = (int) TMath::Min(i, numberIterationsX), r2 = 00900 (int) TMath::Min(i, numberIterationsY); 00901 for (y = r2; y < ssizey - r2; y++) { 00902 for (x = r1; x < ssizex - r1; x++) { 00903 a = spectrum[x][y]; 00904 b = -(spectrum[x - r1][y - r2] + 00905 spectrum[x - r1][y + r2] + spectrum[x + r1][y - 00906 r2] 00907 + spectrum[x + r1][y + r2]) / 4 + 00908 (spectrum[x][y - r2] + spectrum[x - r1][y] + 00909 spectrum[x + r1][y] + spectrum[x][y + r2]) / 2; 00910 if (b < a && b > 0) 00911 a = b; 00912 working_space[x][y] = a; 00913 } 00914 } 00915 for (y = i; y < ssizey - i; y++) { 00916 for (x = i; x < ssizex - i; x++) { 00917 spectrum[x][y] = working_space[x][y]; 00918 } 00919 } 00920 } 00921 } 00922 } 00923 00924 else if (direction == kBackDecreasingWindow) { 00925 if (filterType == kBackSuccessiveFiltering) { 00926 for (i = sampling; i >= 1; i--) { 00927 r1 = (int) TMath::Min(i, numberIterationsX), r2 = 00928 (int) TMath::Min(i, numberIterationsY); 00929 for (y = r2; y < ssizey - r2; y++) { 00930 for (x = r1; x < ssizex - r1; x++) { 00931 a = spectrum[x][y]; 00932 p1 = spectrum[x - r1][y - r2]; 00933 p2 = spectrum[x - r1][y + r2]; 00934 p3 = spectrum[x + r1][y - r2]; 00935 p4 = spectrum[x + r1][y + r2]; 00936 s1 = spectrum[x][y - r2]; 00937 s2 = spectrum[x - r1][y]; 00938 s3 = spectrum[x + r1][y]; 00939 s4 = spectrum[x][y + r2]; 00940 b = (p1 + p2) / 2.0; 00941 if (b > s2) 00942 s2 = b; 00943 b = (p1 + p3) / 2.0; 00944 if (b > s1) 00945 s1 = b; 00946 b = (p2 + p4) / 2.0; 00947 if (b > s4) 00948 s4 = b; 00949 b = (p3 + p4) / 2.0; 00950 if (b > s3) 00951 s3 = b; 00952 s1 = s1 - (p1 + p3) / 2.0; 00953 s2 = s2 - (p1 + p2) / 2.0; 00954 s3 = s3 - (p3 + p4) / 2.0; 00955 s4 = s4 - (p2 + p4) / 2.0; 00956 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + 00957 p3 + 00958 p4) / 4.0; 00959 if (b < a && b > 0) 00960 a = b; 00961 working_space[x][y] = a; 00962 } 00963 } 00964 for (y = r2; y < ssizey - r2; y++) { 00965 for (x = r1; x < ssizex - r1; x++) { 00966 spectrum[x][y] = working_space[x][y]; 00967 } 00968 } 00969 } 00970 } 00971 00972 else if (filterType == kBackOneStepFiltering) { 00973 for (i = sampling; i >= 1; i--) { 00974 r1 = (int) TMath::Min(i, numberIterationsX), r2 = 00975 (int) TMath::Min(i, numberIterationsY); 00976 for (y = r2; y < ssizey - r2; y++) { 00977 for (x = r1; x < ssizex - r1; x++) { 00978 a = spectrum[x][y]; 00979 b = -(spectrum[x - r1][y - r2] + 00980 spectrum[x - r1][y + r2] + spectrum[x + r1][y - 00981 r2] 00982 + spectrum[x + r1][y + r2]) / 4 + 00983 (spectrum[x][y - r2] + spectrum[x - r1][y] + 00984 spectrum[x + r1][y] + spectrum[x][y + r2]) / 2; 00985 if (b < a && b > 0) 00986 a = b; 00987 working_space[x][y] = a; 00988 } 00989 } 00990 for (y = i; y < ssizey - i; y++) { 00991 for (x = i; x < ssizex - i; x++) { 00992 spectrum[x][y] = working_space[x][y]; 00993 } 00994 } 00995 } 00996 } 00997 } 00998 for (i = 0; i < ssizex; i++) 00999 delete[]working_space[i]; 01000 delete[]working_space; 01001 return 0; 01002 } 01003 01004 //_____________________________________________________________________________ 01005 const char* TSpectrum2::SmoothMarkov(float **source, Int_t ssizex, Int_t ssizey, Int_t averWindow) 01006 { 01007 ///////////////////////////////////////////////////////////////////////////// 01008 // TWO-DIMENSIONAL MARKOV SPECTRUM SMOOTHING FUNCTION 01009 // 01010 // This function calculates smoothed spectrum from source spectrum 01011 // based on Markov chain method. 01012 // The result is placed in the array pointed by source pointer. 01013 // 01014 // Function parameters: 01015 // source-pointer to the array of source spectrum 01016 // ssizex-x length of source 01017 // ssizey-y length of source 01018 // averWindow-width of averaging smoothing window 01019 // 01020 ///////////////////////////////////////////////////////////////////////////// 01021 //Begin_Html <!-- 01022 /* --> 01023 <div class=Section1> 01024 01025 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:20.0pt'>Smoothing</span></b></p> 01026 01027 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p> 01028 01029 <p class=MsoNormal><i><span style='font-size:18.0pt'>Goal: Suppression of 01030 statistical fluctuations</span></i></p> 01031 01032 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 01033 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 01034 </span></span><span style='font-size:16.0pt'>the algorithm is based on discrete 01035 Markov chain, which has very simple invariant distribution</span></p> 01036 01037 <p class=MsoNormal> </p> 01038 01039 <p class=MsoNormal> <sub><img width=551 height=63 01040 src="gif/TSpectrum2_Smoothing1.gif"></sub><span style='font-size:16.0pt; 01041 font-family:Arial'> </span></p> 01042 01043 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 01044 font-family:Arial'> </span><sub><img width=28 height=36 01045 src="gif/TSpectrum2_Smoothing2.gif"></sub><span style='font-size:16.0pt; 01046 font-family:Arial'> being defined from the normalization condition </span><sub><img 01047 width=70 height=52 src="gif/TSpectrum2_Smoothing3.gif"></sub></p> 01048 01049 <p class=MsoNormal> </p> 01050 01051 <p class=MsoNormal><span style='font-size:16.0pt;font-family:Arial'> n 01052 is the length of the smoothed spectrum and </span></p> 01053 01054 <p class=MsoNormal> 01055 01056 <table cellpadding=0 cellspacing=0 align=left> 01057 <tr> 01058 <td width=57 height=15></td> 01059 </tr> 01060 <tr> 01061 <td></td> 01062 <td><img width=258 height=60 src="gif/TSpectrum2_Smoothing4.gif"></td> 01063 </tr> 01064 </table> 01065 01066 <span style='font-size:16.0pt;font-family:Arial'> </span></p> 01067 01068 <p class=MsoNormal><span style='font-size:18.0pt'> </span></p> 01069 01070 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 01071 01072 <br clear=ALL> 01073 01074 <p class=MsoNormal style='margin-left:34.2pt;text-align:justify'><span 01075 style='font-size:16.0pt;font-family:Arial'>is the probability of the change of 01076 the peak position from channel i to the channel i+1. </span> <sub><img 01077 width=28 height=36 src="gif/TSpectrum2_Smoothing5.gif"></sub><span 01078 style='font-size:16.0pt;font-family:Arial'>is the normalization constant so 01079 that </span><span style='font-family:Arial'><sub><img width=133 height=34 01080 src="gif/TSpectrum2_Smoothing6.gif"></sub> </span><span style='font-size:16.0pt; 01081 font-family:Arial'>and m is a width of smoothing window. We have extended this 01082 algortihm to two dimensions. </span></p> 01083 01084 <p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p> 01085 01086 <p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p> 01087 01088 <p class=MsoNormal><b><span style='font-size:18.0pt'>const <a 01089 href="http://root.cern.ch/root/html/ListOfTypes.html#char" target="_parent">char</a>* 01090 </span></b><span style='font-size:18.0pt'><a 01091 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b>TSpectrum2::SmoothMarkov</b></a><b>(<a 01092 href="http://root.cern.ch/root/html/ListOfTypes.html#float" target="_parent">float</a> 01093 **fSpectrum, <a href="http://root.cern.ch/root/html/ListOfTypes.html#int" 01094 target="_parent">int</a> ssizex, <a 01095 href="http://root.cern.ch/root/html/ListOfTypes.html#int" target="_parent">int</a> 01096 ssizey, <a href="http://root.cern.ch/root/html/ListOfTypes.html#int" 01097 target="_parent">int</a> averWindow) </b></span></p> 01098 01099 <p class=MsoNormal> </p> 01100 01101 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This 01102 function calculates smoothed spectrum from the source spectrum based on Markov 01103 chain method. The result is placed in the vector pointed by source pointer. On 01104 successful completion it returns 0. On error it returns pointer to the string 01105 describing error.</span></p> 01106 01107 <p class=MsoNormal> </p> 01108 01109 <p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p> 01110 01111 <p class=MsoNormal> <b><span style='font-size:14.0pt'>fSpectrum</span></b>-pointer 01112 to the matrix of source spectrum </p> 01113 01114 <p class=MsoNormal> <b><span style='font-size:14.0pt'>ssizex, ssizey</span></b> 01115 -lengths of the spectrum matrix </p> 01116 01117 <p class=MsoNormal> <b><span style='font-size:14.0pt'>averWindow</span></b>-width 01118 of averaging smoothing window </p> 01119 01120 <p class=MsoNormal> </p> 01121 01122 <p class=MsoNormal><b><i><span style='font-size:18.0pt'>Reference:</span></i></b></p> 01123 01124 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[1] 01125 Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 01126 (1996), 451<b>.</b> </span></p> 01127 01128 </div> 01129 01130 <!-- */ 01131 // --> End_Html 01132 //Begin_Html <!-- 01133 /* --> 01134 <div class=Section1> 01135 01136 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 4 – script Smooth.c 01137 :</span></i></p> 01138 01139 <p class=MsoNormal><span style='font-size:16.0pt'><img width=300 height=209 01140 src="gif/TSpectrum2_Smoothing1.jpg"><img width=297 height=207 01141 src="gif/TSpectrum2_Smoothing2.jpg"></span></p> 01142 01143 <p class=MsoNormal><b><span style='font-size:16.0pt'>Fig. 9 Original noisy 01144 spectrum.</span></b><b><span style='font-size:14.0pt'> </span></b><b><span 01145 style='font-size:16.0pt'>Fig. 10 Smoothed spectrum m=3</span></b></p> 01146 01147 <p class=MsoNormal><b><span style='font-size:16.0pt'>Peaks can hardly be 01148 observed. Peaks become apparent.</span></b></p> 01149 01150 <p class=MsoNormal><b><span style='font-size:14.0pt'><img width=293 height=203 01151 src="gif/TSpectrum2_Smoothing3.jpg"><img width=297 height=205 01152 src="gif/TSpectrum2_Smoothing4.jpg"></span></b></p> 01153 01154 <p class=MsoNormal><b><span style='font-size:16.0pt'>Fig. 11 Smoothed spectrum 01155 m=5 Fig.12 Smoothed spectrum m=7</span></b></p> 01156 01157 <p class=MsoNormal> </p> 01158 01159 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 01160 01161 <p class=MsoNormal>// Example to illustrate the Markov smoothing (class 01162 TSpectrum).</p> 01163 01164 <p class=MsoNormal>// To execute this example, do</p> 01165 01166 <p class=MsoNormal>// root > .x Smooth.C</p> 01167 01168 <p class=MsoNormal>#include <TSpectrum> </p> 01169 01170 <p class=MsoNormal>void Smooth() {</p> 01171 01172 <p class=MsoNormal> Int_t i, j;</p> 01173 01174 <p class=MsoNormal> Double_t nbinsx = 256;</p> 01175 01176 <p class=MsoNormal> Double_t nbinsy = 256; </p> 01177 01178 <p class=MsoNormal> Double_t xmin = 0;</p> 01179 01180 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 01181 01182 <p class=MsoNormal> Double_t ymin = 0;</p> 01183 01184 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 01185 01186 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 01187 01188 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 01189 01190 <p class=MsoNormal> source[i]=new 01191 float[nbinsy]; </p> 01192 01193 <p class=MsoNormal> TH2F *smooth = new 01194 TH2F("smooth","Background 01195 estimation",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 01196 01197 <p class=MsoNormal> TFile *f = new 01198 TFile("spectra2\\TSpectrum2.root");</p> 01199 01200 <p class=MsoNormal> smooth=(TH2F*) f->Get("smooth1;1");</p> 01201 01202 <p class=MsoNormal> TCanvas *Smoothing = new 01203 TCanvas("Smoothing","Markov smoothing",10,10,1000,700);</p> 01204 01205 <p class=MsoNormal> TSpectrum *s = new TSpectrum();</p> 01206 01207 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01208 01209 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 01210 01211 <p class=MsoNormal> source[i][j] = smooth->GetBinContent(i + 01212 1,j + 1); </p> 01213 01214 <p class=MsoNormal> }</p> 01215 01216 <p class=MsoNormal> }</p> 01217 01218 <p class=MsoNormal> s->SmoothMarkov(source,nbinsx,nbinsx,3);//5,7</p> 01219 01220 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01221 01222 <p class=MsoNormal> for (j = 0; j < nbinsy; j++)</p> 01223 01224 <p class=MsoNormal> smooth->SetBinContent(i + 1,j + 1, 01225 source[i][j]); </p> 01226 01227 <p class=MsoNormal> }</p> 01228 01229 <p class=MsoNormal> smooth->Draw("SURF"); </p> 01230 01231 <p class=MsoNormal> }</p> 01232 01233 </div> 01234 01235 <!-- */ 01236 // --> End_Html 01237 int xmin, xmax, ymin, ymax, i, j, l; 01238 double a, b, maxch; 01239 double nom, nip, nim, sp, sm, spx, spy, smx, smy, plocha = 0; 01240 if(averWindow <= 0) 01241 return "Averaging Window must be positive"; 01242 float **working_space = new float* [ssizex]; 01243 for(i = 0; i < ssizex; i++) 01244 working_space[i] = new float[ssizey]; 01245 xmin = 0; 01246 xmax = ssizex - 1; 01247 ymin = 0; 01248 ymax = ssizey - 1; 01249 for(i = 0, maxch = 0; i < ssizex; i++){ 01250 for(j = 0; j < ssizey; j++){ 01251 working_space[i][j] = 0; 01252 if(maxch < source[i][j]) 01253 maxch = source[i][j]; 01254 01255 plocha += source[i][j]; 01256 } 01257 } 01258 if(maxch == 0) { 01259 delete [] working_space; 01260 return 0; 01261 } 01262 01263 nom = 0; 01264 working_space[xmin][ymin] = 1; 01265 for(i = xmin; i < xmax; i++){ 01266 nip = source[i][ymin] / maxch; 01267 nim = source[i + 1][ymin] / maxch; 01268 sp = 0,sm = 0; 01269 for(l = 1; l <= averWindow; l++){ 01270 if((i + l) > xmax) 01271 a = source[xmax][ymin] / maxch; 01272 01273 else 01274 a = source[i + l][ymin] / maxch; 01275 b = a - nip; 01276 if(a + nip <= 0) 01277 a = 1; 01278 01279 else 01280 a = TMath::Sqrt(a + nip); 01281 b = b / a; 01282 b = TMath::Exp(b); 01283 sp = sp + b; 01284 if(i - l + 1 < xmin) 01285 a = source[xmin][ymin] / maxch; 01286 01287 else 01288 a = source[i - l + 1][ymin] / maxch; 01289 b = a - nim; 01290 if(a + nim <= 0) 01291 a = 1; 01292 01293 else 01294 a = TMath::Sqrt(a + nim); 01295 b = b / a; 01296 b = TMath::Exp(b); 01297 sm = sm + b; 01298 } 01299 a = sp / sm; 01300 a = working_space[i + 1][ymin] = a * working_space[i][ymin]; 01301 nom = nom + a; 01302 } 01303 for(i = ymin; i < ymax; i++){ 01304 nip = source[xmin][i] / maxch; 01305 nim = source[xmin][i + 1] / maxch; 01306 sp = 0,sm = 0; 01307 for(l = 1; l <= averWindow; l++){ 01308 if((i + l) > ymax) 01309 a = source[xmin][ymax] / maxch; 01310 01311 else 01312 a = source[xmin][i + l] / maxch; 01313 b = a - nip; 01314 if(a + nip <= 0) 01315 a = 1; 01316 01317 else 01318 a = TMath::Sqrt(a + nip); 01319 b = b / a; 01320 b = TMath::Exp(b); 01321 sp = sp + b; 01322 if(i - l + 1 < ymin) 01323 a = source[xmin][ymin] / maxch; 01324 01325 else 01326 a = source[xmin][i - l + 1] / maxch; 01327 b = a - nim; 01328 if(a + nim <= 0) 01329 a = 1; 01330 01331 else 01332 a = TMath::Sqrt(a + nim); 01333 b = b / a; 01334 b = TMath::Exp(b); 01335 sm = sm + b; 01336 } 01337 a = sp / sm; 01338 a = working_space[xmin][i + 1] = a * working_space[xmin][i]; 01339 nom = nom + a; 01340 } 01341 for(i = xmin; i < xmax; i++){ 01342 for(j = ymin; j < ymax; j++){ 01343 nip = source[i][j + 1] / maxch; 01344 nim = source[i + 1][j + 1] / maxch; 01345 spx = 0,smx = 0; 01346 for(l = 1; l <= averWindow; l++){ 01347 if(i + l > xmax) 01348 a = source[xmax][j] / maxch; 01349 01350 else 01351 a = source[i + l][j] / maxch; 01352 b = a - nip; 01353 if(a + nip <= 0) 01354 a = 1; 01355 01356 else 01357 a = TMath::Sqrt(a + nip); 01358 b = b / a; 01359 b = TMath::Exp(b); 01360 spx = spx + b; 01361 if(i - l + 1 < xmin) 01362 a = source[xmin][j] / maxch; 01363 01364 else 01365 a = source[i - l + 1][j] / maxch; 01366 b = a - nim; 01367 if(a + nim <= 0) 01368 a = 1; 01369 01370 else 01371 a = TMath::Sqrt(a + nim); 01372 b = b / a; 01373 b = TMath::Exp(b); 01374 smx = smx + b; 01375 } 01376 spy = 0,smy = 0; 01377 nip = source[i + 1][j] / maxch; 01378 nim = source[i + 1][j + 1] / maxch; 01379 for (l = 1; l <= averWindow; l++) { 01380 if (j + l > ymax) a = source[i][ymax]/maxch; 01381 else a = source[i][j + l] / maxch; 01382 b = a - nip; 01383 if (a + nip <= 0) a = 1; 01384 else a = TMath::Sqrt(a + nip); 01385 b = b / a; 01386 b = TMath::Exp(b); 01387 spy = spy + b; 01388 if (j - l + 1 < ymin) a = source[i][ymin] / maxch; 01389 else a = source[i][j - l + 1] / maxch; 01390 b = a - nim; 01391 if (a + nim <= 0) a = 1; 01392 else a = TMath::Sqrt(a + nim); 01393 b = b / a; 01394 b = TMath::Exp(b); 01395 smy = smy + b; 01396 } 01397 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx +smy); 01398 working_space[i + 1][j + 1] = a; 01399 nom = nom + a; 01400 } 01401 } 01402 for(i = xmin; i <= xmax; i++){ 01403 for(j = ymin; j <= ymax; j++){ 01404 working_space[i][j] = working_space[i][j] / nom; 01405 } 01406 } 01407 for(i = 0;i < ssizex; i++){ 01408 for(j = 0; j < ssizey; j++){ 01409 source[i][j] = plocha * working_space[i][j]; 01410 } 01411 } 01412 for (i = 0; i < ssizex; i++) 01413 delete[]working_space[i]; 01414 delete[]working_space; 01415 return 0; 01416 } 01417 01418 //______________________________________________________________________________________________________________________________ 01419 const char *TSpectrum2::Deconvolution(float **source, float **resp, 01420 Int_t ssizex, Int_t ssizey, 01421 Int_t numberIterations, 01422 Int_t numberRepetitions, 01423 Double_t boost) 01424 { 01425 ///////////////////////////////////////////////////////////////////////////// 01426 // TWO-DIMENSIONAL DECONVOLUTION FUNCTION 01427 // This function calculates deconvolution from source spectrum 01428 // according to response spectrum 01429 // The result is placed in the matrix pointed by source pointer. 01430 // 01431 // Function parameters: 01432 // source-pointer to the matrix of source spectrum 01433 // resp-pointer to the matrix of response spectrum 01434 // ssizex-x length of source and response spectra 01435 // ssizey-y length of source and response spectra 01436 // numberIterations, for details we refer to manual 01437 // numberRepetitions, for details we refer to manual 01438 // boost, boosting factor, for details we refer to manual 01439 // 01440 ///////////////////////////////////////////////////////////////////////////// 01441 //Begin_Html <!-- 01442 /* --> 01443 <div class=Section1> 01444 01445 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:20.0pt'>Deconvolution</span></b></p> 01446 01447 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p> 01448 01449 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>Goal: 01450 Improvement of the resolution in spectra, decomposition of multiplets</span></i></p> 01451 01452 <p class=MsoNormal><span style='font-size:16.0pt'> </span></p> 01453 01454 <p class=MsoNormal><span style='font-size:16.0pt'>Mathematical formulation of 01455 the 2-dimensional convolution system is</span></p> 01456 01457 <p class=MsoNormal style='margin-left:18.0pt'> 01458 01459 <table cellpadding=0 cellspacing=0 align=left> 01460 <tr> 01461 <td width=0 height=18></td> 01462 </tr> 01463 <tr> 01464 <td></td> 01465 <td><img width=577 height=138 src="gif/TSpectrum2_Deconvolution1.gif"></td> 01466 </tr> 01467 </table> 01468 01469 <span style='font-size:16.0pt'> </span></p> 01470 01471 <p class=MsoNormal><span style='font-size:16.0pt'> </span></p> 01472 01473 <p class=MsoNormal> </p> 01474 01475 <p class=MsoNormal> </p> 01476 01477 <p class=MsoNormal> </p> 01478 01479 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 01480 01481 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 01482 01483 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 01484 01485 <br clear=ALL> 01486 01487 <p class=MsoNormal><span style='font-size:16.0pt'>where h(i,j) is the impulse 01488 response function, x, y are input and output matrices, respectively, <sub><img 01489 width=45 height=24 src="gif/TSpectrum2_Deconvolution2.gif"></sub> are the lengths 01490 of x and h matrices</span><i><span style='font-size:18.0pt'> </span></i></p> 01491 01492 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 01493 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 01494 </span></span><span style='font-size:16.0pt'>let us assume that we know the 01495 response and the output matrices (spectra) of the above given system. </span></p> 01496 01497 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 01498 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 01499 </span></span><span style='font-size:16.0pt'>the deconvolution represents 01500 solution of the overdetermined system of linear equations, i.e., the 01501 calculation of the matrix <b>x.</b></span></p> 01502 01503 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 01504 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 01505 </span></span><span style='font-size:16.0pt'>from numerical stability point of 01506 view the operation of deconvolution is extremely critical (ill-posed problem) 01507 as well as time consuming operation. </span></p> 01508 01509 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 01510 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 01511 </span></span><span style='font-size:16.0pt'>the Gold deconvolution algorithm 01512 proves to work very well even for 2-dimensional systems. Generalization of the 01513 algorithm for 2-dimensional systems was presented in [1], [2].</span></p> 01514 01515 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 01516 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 01517 </span></span><span style='font-size:16.0pt'>for Gold deconvolution algorithm 01518 as well as for boosted deconvolution algorithm we refer also to TSpectrum </span></p> 01519 01520 <p class=MsoNormal><i><span style='font-size:18.0pt'> </span></i></p> 01521 01522 <p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p> 01523 01524 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'>const 01525 <a href="http://root.cern.ch/root/html/ListOfTypes.html#char">char</a>* </span></b><a 01526 name="TSpectrum:Deconvolution1"></a><a 01527 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b><span 01528 style='font-size:18.0pt'>TSpectrum2::Deconvolution</span></b></a><b><span 01529 style='font-size:18.0pt'>(<a 01530 href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> **source, 01531 const <a href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> 01532 **resp, <a href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> ssizex, 01533 <a href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> ssizey, <a 01534 href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> numberIterations, 01535 <a href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> numberRepetitions, 01536 <a href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a></span></b><b><span 01537 style='font-size:16.0pt'> </span></b><b><span style='font-size:18.0pt'>boost)</span></b></p> 01538 01539 <p class=MsoNormal> </p> 01540 01541 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This 01542 function calculates deconvolution from source spectrum according to response 01543 spectrum using Gold deconvolution algorithm. The result is placed in the matrix 01544 pointed by source pointer. On successful completion it returns 0. On error it 01545 returns pointer to the string describing error. If desired after every 01546 numberIterations one can apply boosting operation (exponential function with 01547 exponent given by boost coefficient) and repeat it numberRepetitions times.</span></p> 01548 01549 <p class=MsoNormal> </p> 01550 01551 <p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p> 01552 01553 <p class=MsoNormal> <b><span style='font-size:14.0pt'>source</span></b>-pointer 01554 to the matrix of source spectrum </p> 01555 01556 <p class=MsoNormal> <b><span style='font-size:14.0pt'>resp</span></b>-pointer 01557 to the matrix of response spectrum </p> 01558 01559 <p class=MsoNormal> <b><span style='font-size:14.0pt'>ssizex, ssizey</span></b>-lengths 01560 of the spectrum matrix </p> 01561 01562 <p class=MsoNormal style='text-align:justify'> <b><span 01563 style='font-size:14.0pt'>numberIterations</span></b>-number of iterations </p> 01564 01565 <p class=MsoNormal style='text-align:justify'> <b><span 01566 style='font-size:14.0pt'>numberRepetitions</span></b>-number of repetitions 01567 for boosted deconvolution. It must be </p> 01568 01569 <p class=MsoNormal style='text-align:justify'> greater or equal to one.</p> 01570 01571 <p class=MsoNormal style='text-align:justify'> <b><span 01572 style='font-size:14.0pt'>boost</span></b>-boosting coefficient, applies only 01573 if numberRepetitions is greater than one. </p> 01574 01575 <p class=MsoNormal style='text-align:justify'> <span style='font-size: 01576 14.0pt;color:fuchsia'>Recommended range <1,2>.</span></p> 01577 01578 <p class=MsoNormal> </p> 01579 01580 <p class=MsoNormal><b><i><span style='font-size:18.0pt'>References:</span></i></b></p> 01581 01582 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> [1] 01583 </span><span lang=SK style='font-size:16.0pt'>M. Morháč, J. Kliman, V. 01584 Matoušek, M. Veselský, I. Turzo</span><span style='font-size:16.0pt'>.: 01585 Efficient one- and two-dimensional Gold deconvolution and its application to 01586 gamma-ray spectra decomposition. NIM, A401 (1997) 385-408.</span></p> 01587 01588 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[2] 01589 Morháč M., Matoušek V., Kliman J., Efficient algorithm of multidimensional 01590 deconvolution and its application to nuclear data processing, Digital Signal 01591 Processing 13 (2003) 144. </span></p> 01592 01593 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 01594 01595 </div> 01596 01597 <!-- */ 01598 // --> End_Html 01599 //Begin_Html <!-- 01600 /* --> 01601 <div class=Section1> 01602 01603 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 5 – script Decon.c 01604 :</span></i></p> 01605 01606 <p class=MsoNormal style='margin-left:36.0pt;text-align:justify;text-indent: 01607 -18.0pt'><span style='font-size:16.0pt'>•<span style='font:7.0pt "Times New Roman"'> 01608 </span></span><span style='font-size:16.0pt'>response function (usually peak) 01609 should be shifted to the beginning of the coordinate system (see Fig. 13)</span></p> 01610 01611 <p class=MsoNormal> </p> 01612 01613 <p class=MsoNormal><img width=602 height=418 01614 src="gif/TSpectrum2_Deconvolution1.jpg"></p> 01615 01616 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 01617 13 2-dimensional response spectrum</span></b></p> 01618 01619 <p class=MsoNormal><img width=602 height=418 01620 src="gif/TSpectrum2_Deconvolution2.jpg"></p> 01621 01622 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 01623 14 2-dimensional gamma-gamma-ray input spectrum (before deconvolution)</span></b></p> 01624 01625 <p class=MsoNormal><img width=602 height=418 01626 src="gif/TSpectrum2_Deconvolution3.jpg"></p> 01627 01628 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 01629 15 Spectrum from Fig. 14 after deconvolution (1000 iterations)</span></b></p> 01630 01631 <p class=MsoNormal> </p> 01632 01633 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 01634 01635 <p class=MsoNormal>// Example to illustrate the Gold deconvolution (class 01636 TSpectrum2).</p> 01637 01638 <p class=MsoNormal>// To execute this example, do</p> 01639 01640 <p class=MsoNormal>// root > .x Decon.C</p> 01641 01642 <p class=MsoNormal> </p> 01643 01644 <p class=MsoNormal>#include <TSpectrum2> </p> 01645 01646 <p class=MsoNormal> </p> 01647 01648 <p class=MsoNormal>void Decon() {</p> 01649 01650 <p class=MsoNormal> Int_t i, j;</p> 01651 01652 <p class=MsoNormal> Double_t nbinsx = 256;</p> 01653 01654 <p class=MsoNormal> Double_t nbinsy = 256; </p> 01655 01656 <p class=MsoNormal> Double_t xmin = 0;</p> 01657 01658 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 01659 01660 <p class=MsoNormal> Double_t ymin = 0;</p> 01661 01662 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 01663 01664 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 01665 01666 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 01667 01668 <p class=MsoNormal> source[i]=new 01669 float[nbinsy]; </p> 01670 01671 <p class=MsoNormal> TH2F *decon = new TH2F("decon","Gold 01672 deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 01673 01674 <p class=MsoNormal> TFile *f = new 01675 TFile("spectra2\\TSpectrum2.root");</p> 01676 01677 <p class=MsoNormal> decon=(TH2F*) f->Get("decon1;1");</p> 01678 01679 <p class=MsoNormal> Float_t ** response = new float *[nbinsx]; </p> 01680 01681 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 01682 01683 <p class=MsoNormal> response[i]=new 01684 float[nbinsy]; </p> 01685 01686 <p class=MsoNormal> TH2F *resp = new TH2F("resp","Response 01687 matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 01688 01689 <p class=MsoNormal> resp=(TH2F*) f->Get("resp1;1"); </p> 01690 01691 <p class=MsoNormal> TCanvas *Deconvol = new 01692 TCanvas("Deconvolution","Gold deconvolution",10,10,1000,700);</p> 01693 01694 <p class=MsoNormal> TSpectrum *s = new TSpectrum();</p> 01695 01696 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01697 01698 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 01699 01700 <p class=MsoNormal> source[i][j] = decon->GetBinContent(i + 01701 1,j + 1); </p> 01702 01703 <p class=MsoNormal> }</p> 01704 01705 <p class=MsoNormal> }</p> 01706 01707 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01708 01709 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 01710 01711 <p class=MsoNormal> response[i][j] = resp->GetBinContent(i + 01712 1,j + 1); </p> 01713 01714 <p class=MsoNormal> }</p> 01715 01716 <p class=MsoNormal> } </p> 01717 01718 <p class=MsoNormal> s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1); 01719 </p> 01720 01721 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01722 01723 <p class=MsoNormal> for (j = 0; j < nbinsy; j++)</p> 01724 01725 <p class=MsoNormal> decon->SetBinContent(i + 1,j + 1, source[i][j]); 01726 </p> 01727 01728 <p class=MsoNormal> }</p> 01729 01730 <p class=MsoNormal> </p> 01731 01732 <p class=MsoNormal> decon->Draw("SURF"); </p> 01733 01734 <p class=MsoNormal> }</p> 01735 01736 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 6 – script 01737 Decon2.c :</span></i></p> 01738 01739 <p class=MsoNormal><img width=602 height=418 01740 src="gif/TSpectrum2_Deconvolution4.jpg"></p> 01741 01742 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 01743 16 Response spectrum</span></b></p> 01744 01745 <p class=MsoNormal><img width=602 height=418 01746 src="gif/TSpectrum2_Deconvolution5.jpg"></p> 01747 01748 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 01749 17 Original synthetic input spectrum (before deconvolution). It is composed of 01750 17 peaks. 5 peaks are overlapping in the outlined multiplet and two peaks in 01751 doublet.</span></b></p> 01752 01753 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 01754 width=602 height=418 src="gif/TSpectrum2_Deconvolution6.jpg"></span></b></p> 01755 01756 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 01757 18 Spectrum from Fig. 17 after deconvolution (1000 iterations). Resolution is 01758 improved but the peaks in multiplet remained unresolved.</span></b></p> 01759 01760 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 01761 01762 <p class=MsoNormal>// Example to illustrate the Gold deconvolution (class 01763 TSpectrum2).</p> 01764 01765 <p class=MsoNormal>// To execute this example, do</p> 01766 01767 <p class=MsoNormal>// root > .x Decon2.C</p> 01768 01769 <p class=MsoNormal> </p> 01770 01771 <p class=MsoNormal>#include <TSpectrum2> </p> 01772 01773 <p class=MsoNormal> </p> 01774 01775 <p class=MsoNormal>void Decon2() {</p> 01776 01777 <p class=MsoNormal> Int_t i, j;</p> 01778 01779 <p class=MsoNormal> Double_t nbinsx = 64;</p> 01780 01781 <p class=MsoNormal> Double_t nbinsy = 64; </p> 01782 01783 <p class=MsoNormal> Double_t xmin = 0;</p> 01784 01785 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 01786 01787 <p class=MsoNormal> Double_t ymin = 0;</p> 01788 01789 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 01790 01791 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 01792 01793 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 01794 01795 <p class=MsoNormal> source[i]=new 01796 float[nbinsy]; </p> 01797 01798 <p class=MsoNormal> TH2F *decon = new TH2F("decon","Gold 01799 deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 01800 01801 <p class=MsoNormal> TFile *f = new TFile("spectra2\\TSpectrum2.root");</p> 01802 01803 <p class=MsoNormal> decon=(TH2F*) f->Get("decon2;1");</p> 01804 01805 <p class=MsoNormal> Float_t ** response = new float *[nbinsx]; </p> 01806 01807 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 01808 01809 <p class=MsoNormal> response[i]=new 01810 float[nbinsy]; </p> 01811 01812 <p class=MsoNormal> TH2F *resp = new TH2F("resp","Response 01813 matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 01814 01815 <p class=MsoNormal> resp=(TH2F*) f->Get("resp2;1"); </p> 01816 01817 <p class=MsoNormal> TCanvas *Deconvol = new 01818 TCanvas("Deconvolution","Gold 01819 deconvolution",10,10,1000,700);</p> 01820 01821 <p class=MsoNormal> TSpectrum *s = new TSpectrum();</p> 01822 01823 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01824 01825 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 01826 01827 <p class=MsoNormal> source[i][j] = decon->GetBinContent(i + 01828 1,j + 1); </p> 01829 01830 <p class=MsoNormal> }</p> 01831 01832 <p class=MsoNormal> }</p> 01833 01834 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01835 01836 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 01837 01838 <p class=MsoNormal> response[i][j] = resp->GetBinContent(i + 01839 1,j + 1); </p> 01840 01841 <p class=MsoNormal> }</p> 01842 01843 <p class=MsoNormal> } </p> 01844 01845 <p class=MsoNormal> s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1); 01846 </p> 01847 01848 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01849 01850 <p class=MsoNormal> for (j = 0; j < nbinsy; j++)</p> 01851 01852 <p class=MsoNormal> decon->SetBinContent(i + 1,j + 1, source[i][j]); 01853 </p> 01854 01855 <p class=MsoNormal> }</p> 01856 01857 <p class=MsoNormal> decon->Draw("SURF"); </p> 01858 01859 <p class=MsoNormal> }</p> 01860 01861 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 7 – script 01862 Decon2HR.c :</span></i></p> 01863 01864 <p class=MsoNormal><img width=602 height=418 01865 src="gif/TSpectrum2_Deconvolution7.jpg"></p> 01866 01867 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 01868 19 Spectrum from Fig. 17 after boosted deconvolution (50 iterations repeated 20 01869 times, boosting coefficient was 1.2). All the peaks in multiplet as well as in 01870 doublet are completely decomposed.</span></b></p> 01871 01872 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 01873 01874 <p class=MsoNormal>// Example to illustrate boosted Gold deconvolution (class 01875 TSpectrum2).</p> 01876 01877 <p class=MsoNormal>// To execute this example, do</p> 01878 01879 <p class=MsoNormal>// root > .x Decon2HR.C</p> 01880 01881 <p class=MsoNormal> </p> 01882 01883 <p class=MsoNormal>//#include <TSpectrum2> </p> 01884 01885 <p class=MsoNormal> </p> 01886 01887 <p class=MsoNormal>void Decon2HR() {</p> 01888 01889 <p class=MsoNormal> Int_t i, j;</p> 01890 01891 <p class=MsoNormal> Double_t nbinsx = 64;</p> 01892 01893 <p class=MsoNormal> Double_t nbinsy = 64; </p> 01894 01895 <p class=MsoNormal> Double_t xmin = 0;</p> 01896 01897 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 01898 01899 <p class=MsoNormal> Double_t ymin = 0;</p> 01900 01901 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 01902 01903 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 01904 01905 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 01906 01907 <p class=MsoNormal> source[i]=new 01908 float[nbinsy]; </p> 01909 01910 <p class=MsoNormal> TH2F *decon = new TH2F("decon","Boosted 01911 Gold deconvolution",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 01912 01913 <p class=MsoNormal> TFile *f = new TFile("spectra2\\TSpectrum2.root");</p> 01914 01915 <p class=MsoNormal> decon=(TH2F*) f->Get("decon2;1");</p> 01916 01917 <p class=MsoNormal> Float_t ** response = new float *[nbinsx]; </p> 01918 01919 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 01920 01921 <p class=MsoNormal> response[i]=new 01922 float[nbinsy]; </p> 01923 01924 <p class=MsoNormal> TH2F *resp = new TH2F("resp","Response 01925 matrix",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 01926 01927 <p class=MsoNormal> resp=(TH2F*) f->Get("resp2;1"); </p> 01928 01929 <p class=MsoNormal> TCanvas *Deconvol = new 01930 TCanvas("Deconvolution","Gold 01931 deconvolution",10,10,1000,700);</p> 01932 01933 <p class=MsoNormal> TSpectrum *s = new TSpectrum();</p> 01934 01935 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01936 01937 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 01938 01939 <p class=MsoNormal> source[i][j] = decon->GetBinContent(i + 01940 1,j + 1); </p> 01941 01942 <p class=MsoNormal> }</p> 01943 01944 <p class=MsoNormal> }</p> 01945 01946 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01947 01948 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 01949 01950 <p class=MsoNormal> response[i][j] = resp->GetBinContent(i + 01951 1,j + 1); </p> 01952 01953 <p class=MsoNormal> }</p> 01954 01955 <p class=MsoNormal> } </p> 01956 01957 <p class=MsoNormal> s->Deconvolution(source,response,nbinsx,nbinsy,1000,1,1); 01958 </p> 01959 01960 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 01961 01962 <p class=MsoNormal> for (j = 0; j < nbinsy; j++)</p> 01963 01964 <p class=MsoNormal> decon->SetBinContent(i + 1,j + 1, source[i][j]); 01965 </p> 01966 01967 <p class=MsoNormal> }</p> 01968 01969 <p class=MsoNormal> decon->Draw("SURF"); </p> 01970 01971 <p class=MsoNormal> }</p> 01972 01973 <p class=MsoNormal style='text-align:justify'> </p> 01974 01975 </div> 01976 01977 <!-- */ 01978 // --> End_Html 01979 int i, j, lhx, lhy, i1, i2, j1, j2, k1, k2, lindex, i1min, i1max, 01980 i2min, i2max, j1min, j1max, j2min, j2max, positx = 0, posity = 0, repet; 01981 double lda, ldb, ldc, area, maximum = 0; 01982 if (ssizex <= 0 || ssizey <= 0) 01983 return "Wrong parameters"; 01984 if (numberIterations <= 0) 01985 return "Number of iterations must be positive"; 01986 if (numberRepetitions <= 0) 01987 return "Number of repetitions must be positive"; 01988 double **working_space = new double *[ssizex]; 01989 for (i = 0; i < ssizex; i++) 01990 working_space[i] = new double[5 * ssizey]; 01991 area = 0; 01992 lhx = -1, lhy = -1; 01993 for (i = 0; i < ssizex; i++) { 01994 for (j = 0; j < ssizey; j++) { 01995 lda = resp[i][j]; 01996 if (lda != 0) { 01997 if ((i + 1) > lhx) 01998 lhx = i + 1; 01999 if ((j + 1) > lhy) 02000 lhy = j + 1; 02001 } 02002 working_space[i][j] = lda; 02003 area = area + lda; 02004 if (lda > maximum) { 02005 maximum = lda; 02006 positx = i, posity = j; 02007 } 02008 } 02009 } 02010 if (lhx == -1 || lhy == -1) { 02011 delete [] working_space; 02012 return ("Zero response data"); 02013 } 02014 02015 //calculate ht*y and write into p 02016 for (i2 = 0; i2 < ssizey; i2++) { 02017 for (i1 = 0; i1 < ssizex; i1++) { 02018 ldc = 0; 02019 for (j2 = 0; j2 <= (lhy - 1); j2++) { 02020 for (j1 = 0; j1 <= (lhx - 1); j1++) { 02021 k2 = i2 + j2, k1 = i1 + j1; 02022 if (k2 >= 0 && k2 < ssizey && k1 >= 0 && k1 < ssizex) { 02023 lda = working_space[j1][j2]; 02024 ldb = source[k1][k2]; 02025 ldc = ldc + lda * ldb; 02026 } 02027 } 02028 } 02029 working_space[i1][i2 + ssizey] = ldc; 02030 } 02031 } 02032 02033 //calculate matrix b=ht*h 02034 i1min = -(lhx - 1), i1max = lhx - 1; 02035 i2min = -(lhy - 1), i2max = lhy - 1; 02036 for (i2 = i2min; i2 <= i2max; i2++) { 02037 for (i1 = i1min; i1 <= i1max; i1++) { 02038 ldc = 0; 02039 j2min = -i2; 02040 if (j2min < 0) 02041 j2min = 0; 02042 j2max = lhy - 1 - i2; 02043 if (j2max > lhy - 1) 02044 j2max = lhy - 1; 02045 for (j2 = j2min; j2 <= j2max; j2++) { 02046 j1min = -i1; 02047 if (j1min < 0) 02048 j1min = 0; 02049 j1max = lhx - 1 - i1; 02050 if (j1max > lhx - 1) 02051 j1max = lhx - 1; 02052 for (j1 = j1min; j1 <= j1max; j1++) { 02053 lda = working_space[j1][j2]; 02054 if (i1 + j1 < ssizex && i2 + j2 < ssizey) 02055 ldb = working_space[i1 + j1][i2 + j2]; 02056 else 02057 ldb = 0; 02058 ldc = ldc + lda * ldb; 02059 } 02060 } 02061 working_space[i1 - i1min][i2 - i2min + 2 * ssizey ] = ldc; 02062 } 02063 } 02064 02065 //initialization in x1 matrix 02066 for (i2 = 0; i2 < ssizey; i2++) { 02067 for (i1 = 0; i1 < ssizex; i1++) { 02068 working_space[i1][i2 + 3 * ssizey] = 1; 02069 working_space[i1][i2 + 4 * ssizey] = 0; 02070 } 02071 } 02072 02073 //START OF ITERATIONS 02074 for (repet = 0; repet < numberRepetitions; repet++) { 02075 if (repet != 0) { 02076 for (i = 0; i < ssizex; i++) { 02077 for (j = 0; j < ssizey; j++) { 02078 working_space[i][j + 3 * ssizey] = 02079 TMath::Power(working_space[i][j + 3 * ssizey], boost); 02080 } 02081 } 02082 } 02083 for (lindex = 0; lindex < numberIterations; lindex++) { 02084 for (i2 = 0; i2 < ssizey; i2++) { 02085 for (i1 = 0; i1 < ssizex; i1++) { 02086 ldb = 0; 02087 j2min = i2; 02088 if (j2min > lhy - 1) 02089 j2min = lhy - 1; 02090 j2min = -j2min; 02091 j2max = ssizey - i2 - 1; 02092 if (j2max > lhy - 1) 02093 j2max = lhy - 1; 02094 j1min = i1; 02095 if (j1min > lhx - 1) 02096 j1min = lhx - 1; 02097 j1min = -j1min; 02098 j1max = ssizex - i1 - 1; 02099 if (j1max > lhx - 1) 02100 j1max = lhx - 1; 02101 for (j2 = j2min; j2 <= j2max; j2++) { 02102 for (j1 = j1min; j1 <= j1max; j1++) { 02103 ldc = working_space[j1 - i1min][j2 - i2min + 2 * ssizey]; 02104 lda = working_space[i1 + j1][i2 + j2 + 3 * ssizey]; 02105 ldb = ldb + lda * ldc; 02106 } 02107 } 02108 lda = working_space[i1][i2 + 3 * ssizey]; 02109 ldc = working_space[i1][i2 + 1 * ssizey]; 02110 if (ldc * lda != 0 && ldb != 0) { 02111 lda = lda * ldc / ldb; 02112 } 02113 02114 else 02115 lda = 0; 02116 working_space[i1][i2 + 4 * ssizey] = lda; 02117 } 02118 } 02119 for (i2 = 0; i2 < ssizey; i2++) { 02120 for (i1 = 0; i1 < ssizex; i1++) 02121 working_space[i1][i2 + 3 * ssizey] = 02122 working_space[i1][i2 + 4 * ssizey]; 02123 } 02124 } 02125 } 02126 for (i = 0; i < ssizex; i++) { 02127 for (j = 0; j < ssizey; j++) 02128 source[(i + positx) % ssizex][(j + posity) % ssizey] = 02129 area * working_space[i][j + 3 * ssizey]; 02130 } 02131 for (i = 0; i < ssizex; i++) 02132 delete[]working_space[i]; 02133 delete[]working_space; 02134 return 0; 02135 } 02136 02137 //____________________________________________________________________________ 02138 Int_t TSpectrum2::SearchHighRes(float **source,float **dest, Int_t ssizex, Int_t ssizey, 02139 Double_t sigma, Double_t threshold, 02140 Bool_t backgroundRemove,Int_t deconIterations, 02141 Bool_t markov, Int_t averWindow) 02142 02143 { 02144 ///////////////////////////////////////////////////////////////////////////// 02145 // TWO-DIMENSIONAL HIGH-RESOLUTION PEAK SEARCH FUNCTION // 02146 // This function searches for peaks in source spectrum // 02147 // It is based on deconvolution method. First the background is // 02148 // removed (if desired), then Markov spectrum is calculated // 02149 // (if desired), then the response function is generated // 02150 // according to given sigma and deconvolution is carried out. // 02151 // // 02152 // Function parameters: // 02153 // source-pointer to the matrix of source spectrum // 02154 // dest-pointer to the matrix of resulting deconvolved spectrum // 02155 // ssizex-x length of source spectrum // 02156 // ssizey-y length of source spectrum // 02157 // sigma-sigma of searched peaks, for details we refer to manual // 02158 // threshold-threshold value in % for selected peaks, peaks with // 02159 // amplitude less than threshold*highest_peak/100 // 02160 // are ignored, see manual // 02161 // backgroundRemove-logical variable, set if the removal of // 02162 // background before deconvolution is desired // 02163 // deconIterations-number of iterations in deconvolution operation // 02164 // markov-logical variable, if it is true, first the source spectrum // 02165 // is replaced by new spectrum calculated using Markov // 02166 // chains method. // 02167 // averWindow-averanging window of searched peaks, for details // 02168 // we refer to manual (applies only for Markov method) // 02169 // // 02170 ///////////////////////////////////////////////////////////////////////////// 02171 //Begin_Html <!-- 02172 /* --> 02173 <div class=Section1> 02174 02175 <p class=MsoNormal><b><span style='font-size:20.0pt'>Peaks searching</span></b></p> 02176 02177 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'> </span></i></p> 02178 02179 <p class=MsoNormal style='text-align:justify'><i><span style='font-size:18.0pt'>Goal: 02180 to identify automatically the peaks in spectrum with the presence of the 02181 continuous background, one-fold coincidences (ridges) and statistical 02182 fluctuations - noise.</span></i><span style='font-size:18.0pt'> </span></p> 02183 02184 <p class=MsoNormal><span style='font-size:16.0pt;font-family:Arial'> </span></p> 02185 02186 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 02187 font-family:Arial'>The common problems connected with correct peak 02188 identification in two-dimensional coincidence spectra are</span></p> 02189 02190 <ul style='margin-top:0mm' type=disc> 02191 <li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 02192 font-family:Arial'>non-sensitivity to noise, i.e., only statistically 02193 relevant peaks should be identified</span></li> 02194 <li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 02195 font-family:Arial'>non-sensitivity of the algorithm to continuous 02196 background</span></li> 02197 <li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 02198 font-family:Arial'>non-sensitivity to one-fold coincidences (coincidences 02199 peak – background in both dimensions) and their crossings</span></li> 02200 <li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 02201 font-family:Arial'>ability to identify peaks close to the edges of the 02202 spectrum region. Usually peak finders fail to detect them</span></li> 02203 <li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 02204 font-family:Arial'>resolution, decomposition of doublets and multiplets. 02205 The algorithm should be able to recognize close positioned peaks.</span><span 02206 style='font-size:18.0pt'> </span></li> 02207 <li class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt; 02208 font-family:Arial'>ability to identify peaks with different sigma</span></li> 02209 </ul> 02210 02211 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'> </span></p> 02212 02213 <p class=MsoNormal><i><span style='font-size:18.0pt'>Function:</span></i></p> 02214 02215 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:18.0pt'><a 02216 href="http://root.cern.ch/root/html/ListOfTypes.html#Int_t">Int_t</a> </span></b><a 02217 name="TSpectrum:Search1HighRes"></a><a 02218 href="http://root.cern.ch/root/html/TSpectrum.html#TSpectrum:Fit1Awmi"><b><span 02219 style='font-size:18.0pt'>TSpectrum2::SearchHighRes</span></b></a><b><span 02220 style='font-size:18.0pt'> (<a 02221 href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> **source,<a 02222 href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> **dest, <a 02223 href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> ssizex, <a 02224 href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> ssizey, <a 02225 href="http://root.cern.ch/root/html/ListOfTypes.html#float">float</a> sigma, <a 02226 href="http://root.cern.ch/root/html/ListOfTypes.html#double">double</a> threshold, 02227 <a href="http://root.cern.ch/root/html/ListOfTypes.html#bool">bool</a> backgroundRemove,<a 02228 href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> deconIterations, 02229 <a href="http://root.cern.ch/root/html/ListOfTypes.html#bool">bool</a> markov, 02230 <a href="http://root.cern.ch/root/html/ListOfTypes.html#int">int</a> averWindow) 02231 </span></b></p> 02232 02233 <p class=MsoNormal> </p> 02234 02235 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>This 02236 function searches for peaks in source spectrum. It is based on deconvolution 02237 method. First the background is removed (if desired), then Markov smoothed 02238 spectrum is calculated (if desired), then the response function is generated 02239 according to given sigma and deconvolution is carried out. The order of peaks 02240 is arranged according to their heights in the spectrum after background 02241 elimination. The highest peak is the first in the list. On success it returns 02242 number of found peaks.</span></p> 02243 02244 <p class=MsoNormal> </p> 02245 02246 <p class=MsoNormal><i><span style='font-size:16.0pt;color:red'>Parameters:</span></i></p> 02247 02248 <p class=MsoNormal style='text-align:justify'> <b><span 02249 style='font-size:14.0pt'>source</span></b>-pointer to the matrix of source 02250 spectrum </p> 02251 02252 <p class=MsoNormal style='text-align:justify'> <b><span 02253 style='font-size:14.0pt'>dest</span></b>-resulting spectrum after deconvolution</p> 02254 02255 <p class=MsoNormal style='text-align:justify'> <b><span 02256 style='font-size:14.0pt'>ssizex, ssizey</span></b>-lengths of the source and 02257 destination spectra </p> 02258 02259 <p class=MsoNormal style='text-align:justify'> <b><span 02260 style='font-size:14.0pt'>sigma</span></b>-sigma of searched peaks</p> 02261 02262 <p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span 02263 style='font-size:14.0pt'>threshold</span></b>-<span style='font-size:16.0pt'> </span>threshold 02264 value in % for selected peaks, peaks with amplitude less than 02265 threshold*highest_peak/100 are ignored</p> 02266 02267 <p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span 02268 style='font-size:14.0pt'>backgroundRemove</span></b>-<span style='font-size: 02269 16.0pt'> </span>background_remove-logical variable, true if the removal of 02270 background before deconvolution is desired </p> 02271 02272 <p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span 02273 style='font-size:14.0pt'>deconIterations</span></b>-number of iterations in 02274 deconvolution operation</p> 02275 02276 <p class=MsoNormal style='margin-left:22.8pt;text-align:justify'><b><span 02277 style='font-size:14.0pt'>markov</span></b>-logical variable, if it is true, 02278 first the source spectrum is replaced by new spectrum calculated using Markov 02279 chains method </p> 02280 02281 <p class=MsoNormal style='margin-left:19.95pt;text-align:justify;text-indent: 02282 2.85pt'><b><span style='font-size:14.0pt'>averWindow</span></b>-width of 02283 averaging smoothing window </p> 02284 02285 <p class=MsoNormal> </p> 02286 02287 <p class=MsoNormal><b><i><span style='font-size:18.0pt'>References:</span></i></b></p> 02288 02289 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[1] 02290 M.A. Mariscotti: A method for identification of peaks in the presence of 02291 background and its application to spectrum analysis. NIM 50 (1967), 309-320.</span></p> 02292 02293 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[2] 02294 </span><span lang=SK style='font-size:16.0pt'> M. Morháč, J. Kliman, V. 02295 Matoušek, M. Veselský, I. Turzo</span><span style='font-size:16.0pt'>.:Identification 02296 of peaks in multidimensional coincidence gamma-ray spectra. NIM, A443 (2000) 02297 108-125.</span></p> 02298 02299 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'>[3] 02300 Z.K. Silagadze, A new algorithm for automatic photopeak searches. NIM A 376 02301 (1996), 451.</span></p> 02302 02303 </div> 02304 02305 <!-- */ 02306 // --> End_Html 02307 //Begin_Html <!-- 02308 /* --> 02309 <div class=Section1> 02310 02311 <p class=MsoNormal><b><span style='font-size:18.0pt'>Examples of peak searching 02312 method</span></b></p> 02313 02314 <p class=MsoNormal><span style='font-size:16.0pt'> </span></p> 02315 02316 <p class=MsoNormal style='text-align:justify'><span style='font-size:16.0pt'><a 02317 href="http://root.cern.ch/root/html/src/TSpectrum.cxx.html#TSpectrum:Search1HighRes" 02318 target="_parent">SearchHighRes</a> function provides users with the possibility 02319 to vary the input parameters and with the access to the output deconvolved data 02320 in the destination spectrum. Based on the output data one can tune the 02321 parameters. </span></p> 02322 02323 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 8 – script Src.c:</span></i></p> 02324 02325 <p class=MsoNormal><span style='font-size:16.0pt'><img border=0 width=602 02326 height=455 src="gif/TSpectrum2_Searching1.jpg"></span></p> 02327 02328 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02329 20 Two-dimensional spectrum with found peaks denoted by markers (<sub><img 02330 border=0 width=40 height=19 src="gif/TSpectrum2_Searching2.gif"></sub>, 02331 threshold=5%, 3 iterations steps in the deconvolution)</span></b></p> 02332 02333 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02334 border=0 width=602 height=455 src="gif/TSpectrum2_Searching3.jpg"></span></b></p> 02335 02336 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02337 21 Spectrum from Fig. 20 after background elimination and deconvolution</span></b></p> 02338 02339 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 02340 02341 <p class=MsoNormal>// Example to illustrate high resolution peak searching 02342 function (class TSpectrum).</p> 02343 02344 <p class=MsoNormal>// To execute this example, do</p> 02345 02346 <p class=MsoNormal>// root > .x Src.C</p> 02347 02348 <p class=MsoNormal> </p> 02349 02350 <p class=MsoNormal>#include <TSpectrum2></p> 02351 02352 <p class=MsoNormal> </p> 02353 02354 <p class=MsoNormal>void Src() {</p> 02355 02356 <p class=MsoNormal> Int_t i, j, nfound;</p> 02357 02358 <p class=MsoNormal> Double_t nbinsx = 64;</p> 02359 02360 <p class=MsoNormal> Double_t nbinsy = 64; </p> 02361 02362 <p class=MsoNormal> Double_t xmin = 0;</p> 02363 02364 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 02365 02366 <p class=MsoNormal> Double_t ymin = 0;</p> 02367 02368 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 02369 02370 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 02371 02372 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02373 02374 <p class=MsoNormal> source[i]=new 02375 float[nbinsy];</p> 02376 02377 <p class=MsoNormal> Float_t ** dest = new float *[nbinsx]; </p> 02378 02379 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02380 02381 <p class=MsoNormal> dest[i]=new 02382 float[nbinsy];</p> 02383 02384 <p class=MsoNormal> TH2F *search = new TH2F("search","High 02385 resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 02386 02387 <p class=MsoNormal> TFile *f = new TFile("spectra2\\TSpectrum2.root");</p> 02388 02389 <p class=MsoNormal> search=(TH2F*) f->Get("search4;1");</p> 02390 02391 <p class=MsoNormal> TCanvas *Searching = new 02392 TCanvas("Searching","High resolution peak 02393 searching",10,10,1000,700);</p> 02394 02395 <p class=MsoNormal> TSpectrum2 *s = new TSpectrum2();</p> 02396 02397 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 02398 02399 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 02400 02401 <p class=MsoNormal> source[i][j] = search->GetBinContent(i + 02402 1,j + 1); </p> 02403 02404 <p class=MsoNormal> }</p> 02405 02406 <p class=MsoNormal> } </p> 02407 02408 <p class=MsoNormal> nfound = s->SearchHighRes(source, dest, nbinsx, 02409 nbinsy, 2, 5, kTRUE, 3, kFALSE, 3); </p> 02410 02411 <p class=MsoNormal> printf("Found %d candidate peaks\n",nfound);</p> 02412 02413 <p class=MsoNormal> for(i=0;i<nfound;i++)</p> 02414 02415 <p class=MsoNormal> printf("posx= %d, posy= %d, value= 02416 %d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5), 02417 (int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]); </p> 02418 02419 <p class=MsoNormal>}</p> 02420 02421 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 9 – script Src2.c:</span></i></p> 02422 02423 <p class=MsoNormal><img border=0 width=602 height=455 02424 src="gif/TSpectrum2_Searching4.jpg"></p> 02425 02426 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02427 22 Two-dimensional noisy spectrum with found peaks denoted by markers (<sub><img 02428 border=0 width=40 height=19 src="gif/TSpectrum2_Searching2.gif"></sub>, 02429 threshold=10%, 10 iterations steps in the deconvolution). One can observe that 02430 the algorithm is insensitive to the crossings of one-dimensional ridges. It 02431 identifies only two-cooincidence peaks.</span></b></p> 02432 02433 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02434 border=0 width=602 height=455 src="gif/TSpectrum2_Searching5.jpg"></span></b></p> 02435 02436 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02437 23 Spectrum from Fig. 22 after background elimination and deconvolution</span></b></p> 02438 02439 <p class=MsoNormal> </p> 02440 02441 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 02442 02443 <p class=MsoNormal>// Example to illustrate high resolution peak searching 02444 function (class TSpectrum).</p> 02445 02446 <p class=MsoNormal>// To execute this example, do</p> 02447 02448 <p class=MsoNormal>// root > .x Src2.C</p> 02449 02450 <p class=MsoNormal> </p> 02451 02452 <p class=MsoNormal>#include <TSpectrum2></p> 02453 02454 <p class=MsoNormal> </p> 02455 02456 <p class=MsoNormal>void Src2() {</p> 02457 02458 <p class=MsoNormal> Int_t i, j, nfound;</p> 02459 02460 <p class=MsoNormal> Double_t nbinsx = 256;</p> 02461 02462 <p class=MsoNormal> Double_t nbinsy = 256; </p> 02463 02464 <p class=MsoNormal> Double_t xmin = 0;</p> 02465 02466 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 02467 02468 <p class=MsoNormal> Double_t ymin = 0;</p> 02469 02470 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 02471 02472 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 02473 02474 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02475 02476 <p class=MsoNormal> source[i]=new 02477 float[nbinsy];</p> 02478 02479 <p class=MsoNormal> Float_t ** dest = new float *[nbinsx]; </p> 02480 02481 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02482 02483 <p class=MsoNormal> dest[i]=new 02484 float[nbinsy];</p> 02485 02486 <p class=MsoNormal> TH2F *search = new TH2F("search","High 02487 resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 02488 02489 <p class=MsoNormal> TFile *f = new 02490 TFile("spectra2\\TSpectrum2.root");</p> 02491 02492 <p class=MsoNormal> search=(TH2F*) f->Get("back3;1");</p> 02493 02494 <p class=MsoNormal> TCanvas *Searching = new 02495 TCanvas("Searching","High resolution peak 02496 searching",10,10,1000,700);</p> 02497 02498 <p class=MsoNormal> TSpectrum2 *s = new TSpectrum2();</p> 02499 02500 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 02501 02502 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 02503 02504 <p class=MsoNormal> source[i][j] = search->GetBinContent(i + 02505 1,j + 1); </p> 02506 02507 <p class=MsoNormal> }</p> 02508 02509 <p class=MsoNormal> } </p> 02510 02511 <p class=MsoNormal> nfound = s->SearchHighRes(source, dest, nbinsx, 02512 nbinsy, 2, 10, kTRUE, 10, kFALSE, 3); </p> 02513 02514 <p class=MsoNormal> printf("Found %d candidate peaks\n",nfound);</p> 02515 02516 <p class=MsoNormal> for(i=0;i<nfound;i++)</p> 02517 02518 <p class=MsoNormal> printf("posx= %d, posy= %d, value= 02519 %d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5), 02520 (int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]); </p> 02521 02522 <p class=MsoNormal>}</p> 02523 02524 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 10 – script Src3.c:</span></i></p> 02525 02526 <p class=MsoNormal><img border=0 width=602 height=455 02527 src="gif/TSpectrum2_Searching6.jpg"></p> 02528 02529 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02530 24 Two-dimensional spectrum with 15 found peaks denoted by markers. Some peaks 02531 are positioned close to each other. It is necessary to increase number of 02532 iterations in the deconvolution. In next 3 Figs. we shall study the influence 02533 of this parameter.</span></b></p> 02534 02535 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02536 border=0 width=602 height=455 src="gif/TSpectrum2_Searching7.jpg"></span></b></p> 02537 02538 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02539 25 Spectrum from Fig. 24 after deconvolution (# of iterations = 3). Number of 02540 identified peaks = 13.</span></b></p> 02541 02542 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02543 border=0 width=602 height=455 src="gif/TSpectrum2_Searching8.jpg"></span></b></p> 02544 02545 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02546 26 Spectrum from Fig. 24 after deconvolution (# of iterations = 10). Number of 02547 identified peaks = 13.</span></b></p> 02548 02549 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02550 border=0 width=602 height=455 src="gif/TSpectrum2_Searching9.jpg"></span></b></p> 02551 02552 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02553 27 Spectrum from Fig. 24 after deconvolution (# of iterations = 100). Number of 02554 identified peaks = 15. Now the algorithm is able to decompose two doublets in 02555 the spectrum.</span></b></p> 02556 02557 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'> </span></b></p> 02558 02559 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 02560 02561 <p class=MsoNormal>// Example to illustrate high resolution peak searching 02562 function (class TSpectrum).</p> 02563 02564 <p class=MsoNormal>// To execute this example, do</p> 02565 02566 <p class=MsoNormal>// root > .x Src3.C</p> 02567 02568 <p class=MsoNormal> </p> 02569 02570 <p class=MsoNormal>#include <TSpectrum2></p> 02571 02572 <p class=MsoNormal> </p> 02573 02574 <p class=MsoNormal>void Src3() {</p> 02575 02576 <p class=MsoNormal> Int_t i, j, nfound;</p> 02577 02578 <p class=MsoNormal> Double_t nbinsx = 64;</p> 02579 02580 <p class=MsoNormal> Double_t nbinsy = 64; </p> 02581 02582 <p class=MsoNormal> Double_t xmin = 0;</p> 02583 02584 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 02585 02586 <p class=MsoNormal> Double_t ymin = 0;</p> 02587 02588 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 02589 02590 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 02591 02592 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02593 02594 <p class=MsoNormal> source[i]=new 02595 float[nbinsy];</p> 02596 02597 <p class=MsoNormal> Float_t ** dest = new float *[nbinsx]; </p> 02598 02599 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02600 02601 <p class=MsoNormal> dest[i]=new 02602 float[nbinsy];</p> 02603 02604 <p class=MsoNormal> TH2F *search = new TH2F("search","High 02605 resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 02606 02607 <p class=MsoNormal> TFile *f = new 02608 TFile("spectra2\\TSpectrum2.root");</p> 02609 02610 <p class=MsoNormal> search=(TH2F*) f->Get("search1;1");</p> 02611 02612 <p class=MsoNormal> TCanvas *Searching = new 02613 TCanvas("Searching","High resolution peak 02614 searching",10,10,1000,700);</p> 02615 02616 <p class=MsoNormal> TSpectrum2 *s = new TSpectrum2();</p> 02617 02618 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 02619 02620 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 02621 02622 <p class=MsoNormal> source[i][j] = search->GetBinContent(i + 02623 1,j + 1); </p> 02624 02625 <p class=MsoNormal> }</p> 02626 02627 <p class=MsoNormal> } </p> 02628 02629 <p class=MsoNormal> nfound = s->SearchHighRes(source, dest, nbinsx, 02630 nbinsy, 2, 2, kFALSE, 3, kFALSE, 1);//3, 10, 100 </p> 02631 02632 <p class=MsoNormal> printf("Found %d candidate peaks\n",nfound);</p> 02633 02634 <p class=MsoNormal> </p> 02635 02636 <p class=MsoNormal> for(i=0;i<nfound;i++)</p> 02637 02638 <p class=MsoNormal> printf("posx= %d, posy= %d, value= 02639 %d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5), 02640 (int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]); </p> 02641 02642 <p class=MsoNormal>}</p> 02643 02644 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 11 – script Src4.c:</span></i></p> 02645 02646 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02647 border=0 width=602 height=455 src="gif/TSpectrum2_Searching10.jpg"></span></b></p> 02648 02649 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02650 28 Two-dimensional spectrum with peaks with different sigma denoted by markers (<sub><img 02651 border=0 width=39 height=19 src="gif/TSpectrum2_Searching11.gif"></sub>, 02652 threshold=5%, 10 iterations steps in the deconvolution, Markov smoothing with 02653 window=3)</span></b></p> 02654 02655 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02656 border=0 width=602 height=455 src="gif/TSpectrum2_Searching12.jpg"></span></b></p> 02657 02658 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02659 29 Spectrum from Fig. 28 after smoothing and deconvolution.</span></b></p> 02660 02661 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'> </span></b></p> 02662 02663 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 02664 02665 <p class=MsoNormal>// Example to illustrate high resolution peak searching 02666 function (class TSpectrum).</p> 02667 02668 <p class=MsoNormal>// To execute this example, do</p> 02669 02670 <p class=MsoNormal>// root > .x Src4.C</p> 02671 02672 <p class=MsoNormal> </p> 02673 02674 <p class=MsoNormal>#include <TSpectrum2></p> 02675 02676 <p class=MsoNormal> </p> 02677 02678 <p class=MsoNormal>void Src4() {</p> 02679 02680 <p class=MsoNormal> Int_t i, j, nfound;</p> 02681 02682 <p class=MsoNormal> Double_t nbinsx = 64;</p> 02683 02684 <p class=MsoNormal> Double_t nbinsy = 64; </p> 02685 02686 <p class=MsoNormal> Double_t xmin = 0;</p> 02687 02688 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 02689 02690 <p class=MsoNormal> Double_t ymin = 0;</p> 02691 02692 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 02693 02694 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 02695 02696 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02697 02698 <p class=MsoNormal> source[i]=new 02699 float[nbinsy];</p> 02700 02701 <p class=MsoNormal> Float_t ** dest = new float *[nbinsx]; </p> 02702 02703 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02704 02705 <p class=MsoNormal> dest[i]=new 02706 float[nbinsy];</p> 02707 02708 <p class=MsoNormal> TH2F *search = new TH2F("search","High 02709 resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 02710 02711 <p class=MsoNormal> TFile *f = new 02712 TFile("spectra2\\TSpectrum2.root");</p> 02713 02714 <p class=MsoNormal> search=(TH2F*) f->Get("search2;1");</p> 02715 02716 <p class=MsoNormal> TCanvas *Searching = new 02717 TCanvas("Searching","High resolution peak 02718 searching",10,10,1000,700);</p> 02719 02720 <p class=MsoNormal> TSpectrum2 *s = new TSpectrum2();</p> 02721 02722 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 02723 02724 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 02725 02726 <p class=MsoNormal> source[i][j] = search->GetBinContent(i + 02727 1,j + 1); </p> 02728 02729 <p class=MsoNormal> }</p> 02730 02731 <p class=MsoNormal> } </p> 02732 02733 <p class=MsoNormal> nfound = s->SearchHighRes(source, dest, nbinsx, 02734 nbinsy, 3, 5, kFALSE, 10, kTRUE, 3); </p> 02735 02736 <p class=MsoNormal> printf("Found %d candidate peaks\n",nfound);</p> 02737 02738 <p class=MsoNormal> for(i=0;i<nfound;i++)</p> 02739 02740 <p class=MsoNormal> printf("posx= %d, posy= %d, value= 02741 %d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5), 02742 (int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]); </p> 02743 02744 <p class=MsoNormal>}</p> 02745 02746 <p class=MsoNormal><i><span style='font-size:16.0pt'>Example 12 – script Src5.c:</span></i></p> 02747 02748 <p class=MsoNormal><img border=0 width=602 height=455 02749 src="gif/TSpectrum2_Searching13.jpg"></p> 02750 02751 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02752 30 Two-dimensional spectrum with peaks positioned close to the edges denoted by 02753 markers (<sub><img border=0 width=40 height=19 02754 src="gif/TSpectrum2_Searching2.gif"></sub>, threshold=5%, 10 iterations 02755 steps in the deconvolution)</span></b></p> 02756 02757 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'><img 02758 border=0 width=602 height=455 src="gif/TSpectrum2_Searching14.jpg"></span></b></p> 02759 02760 <p class=MsoNormal style='text-align:justify'><b><span style='font-size:16.0pt'>Fig. 02761 31 Spectrum from Fig. 30 after deconvolution.</span></b></p> 02762 02763 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'> </span></b></p> 02764 02765 <p class=MsoNormal><b><span style='font-size:16.0pt;color:#339966'>Script:</span></b></p> 02766 02767 <p class=MsoNormal>// Example to illustrate high resolution peak searching 02768 function (class TSpectrum).</p> 02769 02770 <p class=MsoNormal>// To execute this example, do</p> 02771 02772 <p class=MsoNormal>// root > .x Src5.C</p> 02773 02774 <p class=MsoNormal> </p> 02775 02776 <p class=MsoNormal>#include <TSpectrum2></p> 02777 02778 <p class=MsoNormal> </p> 02779 02780 <p class=MsoNormal>void Src5() {</p> 02781 02782 <p class=MsoNormal> Int_t i, j, nfound;</p> 02783 02784 <p class=MsoNormal> Double_t nbinsx = 64;</p> 02785 02786 <p class=MsoNormal> Double_t nbinsy = 64; </p> 02787 02788 <p class=MsoNormal> Double_t xmin = 0;</p> 02789 02790 <p class=MsoNormal> Double_t xmax = (Double_t)nbinsx;</p> 02791 02792 <p class=MsoNormal> Double_t ymin = 0;</p> 02793 02794 <p class=MsoNormal> Double_t ymax = (Double_t)nbinsy; </p> 02795 02796 <p class=MsoNormal> Float_t ** source = new float *[nbinsx]; </p> 02797 02798 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02799 02800 <p class=MsoNormal> source[i]=new 02801 float[nbinsy];</p> 02802 02803 <p class=MsoNormal> Float_t ** dest = new float *[nbinsx]; </p> 02804 02805 <p class=MsoNormal> for (i=0;i<nbinsx;i++)</p> 02806 02807 <p class=MsoNormal> dest[i]=new 02808 float[nbinsy];</p> 02809 02810 <p class=MsoNormal> TH2F *search = new TH2F("search","High 02811 resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);</p> 02812 02813 <p class=MsoNormal> TFile *f = new 02814 TFile("spectra2\\TSpectrum2.root");</p> 02815 02816 <p class=MsoNormal> search=(TH2F*) f->Get("search3;1");</p> 02817 02818 <p class=MsoNormal> TCanvas *Searching = new 02819 TCanvas("Searching","High resolution peak 02820 searching",10,10,1000,700);</p> 02821 02822 <p class=MsoNormal> TSpectrum2 *s = new TSpectrum2();</p> 02823 02824 <p class=MsoNormal> for (i = 0; i < nbinsx; i++){</p> 02825 02826 <p class=MsoNormal> for (j = 0; j < nbinsy; j++){</p> 02827 02828 <p class=MsoNormal> source[i][j] = search->GetBinContent(i + 02829 1,j + 1); </p> 02830 02831 <p class=MsoNormal> }</p> 02832 02833 <p class=MsoNormal> } </p> 02834 02835 <p class=MsoNormal> nfound = s->SearchHighRes(source, dest, nbinsx, 02836 nbinsy, 2, 5, kFALSE, 10, kFALSE, 1); </p> 02837 02838 <p class=MsoNormal> printf("Found %d candidate peaks\n",nfound);</p> 02839 02840 <p class=MsoNormal> for(i=0;i<nfound;i++)</p> 02841 02842 <p class=MsoNormal> printf("posx= %d, posy= %d, value= 02843 %d\n",(int)(fPositionX[i]+0.5), (int)(fPositionY[i]+0.5), 02844 (int)source[(int)(fPositionX[i]+0.5)][(int)(fPositionY[i]+0.5)]); </p> 02845 02846 <p class=MsoNormal>}</p> 02847 02848 </div> 02849 02850 <!-- */ 02851 // --> End_Html 02852 int number_of_iterations = (int)(4 * sigma + 0.5); 02853 int k, lindex, priz; 02854 double lda, ldb, ldc, area, maximum; 02855 int xmin, xmax, l, peak_index = 0, ssizex_ext = ssizex + 4 * number_of_iterations, ssizey_ext = ssizey + 4 * number_of_iterations, shift = 2 * number_of_iterations; 02856 int ymin, ymax, i, j; 02857 double a, b, ax, ay, maxch, plocha = 0; 02858 double nom, nip, nim, sp, sm, spx, spy, smx, smy; 02859 double p1, p2, p3, p4, s1, s2, s3, s4; 02860 int x, y; 02861 int lhx, lhy, i1, i2, j1, j2, k1, k2, i1min, i1max, i2min, i2max, j1min, j1max, j2min, j2max, positx, posity; 02862 if (sigma < 1) { 02863 Error("SearchHighRes", "Invalid sigma, must be greater than or equal to 1"); 02864 return 0; 02865 } 02866 02867 if(threshold<=0||threshold>=100){ 02868 Error("SearchHighRes", "Invalid threshold, must be positive and less than 100"); 02869 return 0; 02870 } 02871 02872 j = (int) (5.0 * sigma + 0.5); 02873 if (j >= PEAK_WINDOW / 2) { 02874 Error("SearchHighRes", "Too large sigma"); 02875 return 0; 02876 } 02877 02878 if (markov == true) { 02879 if (averWindow <= 0) { 02880 Error("SearchHighRes", "Averanging window must be positive"); 02881 return 0; 02882 } 02883 } 02884 if(backgroundRemove == true){ 02885 if(ssizex_ext < 2 * number_of_iterations + 1 || ssizey_ext < 2 * number_of_iterations + 1){ 02886 Error("SearchHighRes", "Too large clipping window"); 02887 return 0; 02888 } 02889 } 02890 i = (int)(4 * sigma + 0.5); 02891 i = 4 * i; 02892 double **working_space = new double *[ssizex + i]; 02893 for (j = 0; j < ssizex + i; j++) { 02894 Double_t *wsk = working_space[j] = new double[16 * (ssizey + i)]; 02895 for (k=0;k<16 * (ssizey + i);k++) wsk[k] = 0; 02896 } 02897 for(j = 0; j < ssizey_ext; j++){ 02898 for(i = 0; i < ssizex_ext; i++){ 02899 if(i < shift){ 02900 if(j < shift) 02901 working_space[i][j + ssizey_ext] = source[0][0]; 02902 02903 else if(j >= ssizey + shift) 02904 working_space[i][j + ssizey_ext] = source[0][ssizey - 1]; 02905 02906 else 02907 working_space[i][j + ssizey_ext] = source[0][j - shift]; 02908 } 02909 02910 else if(i >= ssizex + shift){ 02911 if(j < shift) 02912 working_space[i][j + ssizey_ext] = source[ssizex - 1][0]; 02913 02914 else if(j >= ssizey + shift) 02915 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1]; 02916 02917 else 02918 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift]; 02919 } 02920 02921 else{ 02922 if(j < shift) 02923 working_space[i][j + ssizey_ext] = source[i - shift][0]; 02924 02925 else if(j >= ssizey + shift) 02926 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1]; 02927 02928 else 02929 working_space[i][j + ssizey_ext] = source[i - shift][j - shift]; 02930 } 02931 } 02932 } 02933 if(backgroundRemove == true){ 02934 for(i = 1; i <= number_of_iterations; i++){ 02935 for(y = i; y < ssizey_ext - i; y++){ 02936 for(x = i; x < ssizex_ext - i; x++){ 02937 a = working_space[x][y + ssizey_ext]; 02938 p1 = working_space[x - i][y + ssizey_ext - i]; 02939 p2 = working_space[x - i][y + ssizey_ext + i]; 02940 p3 = working_space[x + i][y + ssizey_ext - i]; 02941 p4 = working_space[x + i][y + ssizey_ext + i]; 02942 s1 = working_space[x][y + ssizey_ext - i]; 02943 s2 = working_space[x - i][y + ssizey_ext]; 02944 s3 = working_space[x + i][y + ssizey_ext]; 02945 s4 = working_space[x][y + ssizey_ext + i]; 02946 b = (p1 + p2) / 2.0; 02947 if(b > s2) 02948 s2 = b; 02949 b = (p1 + p3) / 2.0; 02950 if(b > s1) 02951 s1 = b; 02952 b = (p2 + p4) / 2.0; 02953 if(b > s4) 02954 s4 = b; 02955 b = (p3 + p4) / 2.0; 02956 if(b > s3) 02957 s3 = b; 02958 s1 = s1 - (p1 + p3) / 2.0; 02959 s2 = s2 - (p1 + p2) / 2.0; 02960 s3 = s3 - (p3 + p4) / 2.0; 02961 s4 = s4 - (p2 + p4) / 2.0; 02962 b = (s1 + s4) / 2.0 + (s2 + s3) / 2.0 + (p1 + p2 + p3 + p4) / 4.0; 02963 if(b < a) 02964 a = b; 02965 working_space[x][y] = a; 02966 } 02967 } 02968 for(y = i;y < ssizey_ext - i; y++){ 02969 for(x = i; x < ssizex_ext - i; x++){ 02970 working_space[x][y + ssizey_ext] = working_space[x][y]; 02971 } 02972 } 02973 } 02974 for(j = 0;j < ssizey_ext; j++){ 02975 for(i = 0; i < ssizex_ext; i++){ 02976 if(i < shift){ 02977 if(j < shift) 02978 working_space[i][j + ssizey_ext] = source[0][0] - working_space[i][j + ssizey_ext]; 02979 02980 else if(j >= ssizey + shift) 02981 working_space[i][j + ssizey_ext] = source[0][ssizey - 1] - working_space[i][j + ssizey_ext]; 02982 02983 else 02984 working_space[i][j + ssizey_ext] = source[0][j - shift] - working_space[i][j + ssizey_ext]; 02985 } 02986 02987 else if(i >= ssizex + shift){ 02988 if(j < shift) 02989 working_space[i][j + ssizey_ext] = source[ssizex - 1][0] - working_space[i][j + ssizey_ext]; 02990 02991 else if(j >= ssizey + shift) 02992 working_space[i][j + ssizey_ext] = source[ssizex - 1][ssizey - 1] - working_space[i][j + ssizey_ext]; 02993 02994 else 02995 working_space[i][j + ssizey_ext] = source[ssizex - 1][j - shift] - working_space[i][j + ssizey_ext]; 02996 } 02997 02998 else{ 02999 if(j < shift) 03000 working_space[i][j + ssizey_ext] = source[i - shift][0] - working_space[i][j + ssizey_ext]; 03001 03002 else if(j >= ssizey + shift) 03003 working_space[i][j + ssizey_ext] = source[i - shift][ssizey - 1] - working_space[i][j + ssizey_ext]; 03004 03005 else 03006 working_space[i][j + ssizey_ext] = source[i - shift][j - shift] - working_space[i][j + ssizey_ext]; 03007 } 03008 } 03009 } 03010 } 03011 for(j = 0; j < ssizey_ext; j++){ 03012 for(i = 0; i < ssizex_ext; i++){ 03013 working_space[i][j + 15*ssizey_ext] = working_space[i][j + ssizey_ext]; 03014 } 03015 } 03016 if(markov == true){ 03017 for(i = 0;i < ssizex_ext; i++){ 03018 for(j = 0; j < ssizey_ext; j++) 03019 working_space[i][j + 2 * ssizex_ext] = working_space[i][ssizey_ext + j]; 03020 } 03021 xmin = 0; 03022 xmax = ssizex_ext - 1; 03023 ymin = 0; 03024 ymax = ssizey_ext - 1; 03025 for(i = 0, maxch = 0; i < ssizex_ext; i++){ 03026 for(j = 0; j < ssizey_ext; j++){ 03027 working_space[i][j] = 0; 03028 if(maxch < working_space[i][j + 2 * ssizey_ext]) 03029 maxch = working_space[i][j + 2 * ssizey_ext]; 03030 plocha += working_space[i][j + 2 * ssizey_ext]; 03031 } 03032 } 03033 if(maxch == 0) { 03034 delete [] working_space; 03035 return 0; 03036 } 03037 03038 nom=0; 03039 working_space[xmin][ymin] = 1; 03040 for(i = xmin; i < xmax; i++){ 03041 nip = working_space[i][ymin + 2 * ssizey_ext] / maxch; 03042 nim = working_space[i + 1][ymin + 2 * ssizey_ext] / maxch; 03043 sp = 0,sm = 0; 03044 for(l = 1;l <= averWindow; l++){ 03045 if((i + l) > xmax) 03046 a = working_space[xmax][ymin + 2 * ssizey_ext] / maxch; 03047 03048 else 03049 a = working_space[i + l][ymin + 2 * ssizey_ext] / maxch; 03050 03051 b = a - nip; 03052 if(a + nip <= 0) 03053 a = 1; 03054 03055 else 03056 a=TMath::Sqrt(a + nip); 03057 b = b / a; 03058 b = TMath::Exp(b); 03059 sp = sp + b; 03060 if(i - l + 1 < xmin) 03061 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch; 03062 03063 else 03064 a = working_space[i - l + 1][ymin + 2 * ssizey_ext] / maxch; 03065 b = a - nim; 03066 if(a + nim <= 0) 03067 a = 1; 03068 03069 else 03070 a=TMath::Sqrt(a + nim); 03071 b = b / a; 03072 b = TMath::Exp(b); 03073 sm = sm + b; 03074 } 03075 a = sp / sm; 03076 a = working_space[i + 1][ymin] = a * working_space[i][ymin]; 03077 nom = nom + a; 03078 } 03079 for(i = ymin; i < ymax; i++){ 03080 nip = working_space[xmin][i + 2 * ssizey_ext] / maxch; 03081 nim = working_space[xmin][i + 1 + 2 * ssizey_ext] / maxch; 03082 sp = 0,sm = 0; 03083 for(l = 1; l <= averWindow; l++){ 03084 if((i + l) > ymax) 03085 a = working_space[xmin][ymax + 2 * ssizey_ext] / maxch; 03086 03087 else 03088 a = working_space[xmin][i + l + 2 * ssizey_ext] / maxch; 03089 b = a - nip; 03090 if(a + nip <= 0) 03091 a=1; 03092 03093 else 03094 a=TMath::Sqrt(a + nip); 03095 b = b / a; 03096 b = TMath::Exp(b); 03097 sp = sp + b; 03098 if(i - l + 1 < ymin) 03099 a = working_space[xmin][ymin + 2 * ssizey_ext] / maxch; 03100 03101 else 03102 a = working_space[xmin][i - l + 1 + 2 * ssizey_ext] / maxch; 03103 b = a - nim; 03104 if(a + nim <= 0) 03105 a = 1; 03106 03107 else 03108 a=TMath::Sqrt(a + nim); 03109 b = b / a; 03110 b = TMath::Exp(b); 03111 sm = sm + b; 03112 } 03113 a = sp / sm; 03114 a = working_space[xmin][i + 1] = a * working_space[xmin][i]; 03115 nom = nom + a; 03116 } 03117 for(i = xmin; i < xmax; i++){ 03118 for(j = ymin; j < ymax; j++){ 03119 nip = working_space[i][j + 1 + 2 * ssizey_ext] / maxch; 03120 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch; 03121 spx = 0,smx = 0; 03122 for(l = 1; l <= averWindow; l++){ 03123 if(i + l > xmax) 03124 a = working_space[xmax][j + 2 * ssizey_ext] / maxch; 03125 03126 else 03127 a = working_space[i + l][j + 2 * ssizey_ext] / maxch; 03128 b = a - nip; 03129 if(a + nip <= 0) 03130 a = 1; 03131 03132 else 03133 a=TMath::Sqrt(a + nip); 03134 b = b / a; 03135 b = TMath::Exp(b); 03136 spx = spx + b; 03137 if(i - l + 1 < xmin) 03138 a = working_space[xmin][j + 2 * ssizey_ext] / maxch; 03139 03140 else 03141 a = working_space[i - l + 1][j + 2 * ssizey_ext] / maxch; 03142 b = a - nim; 03143 if(a + nim <= 0) 03144 a=1; 03145 03146 else 03147 a=TMath::Sqrt(a + nim); 03148 b = b / a; 03149 b = TMath::Exp(b); 03150 smx = smx + b; 03151 } 03152 spy = 0,smy = 0; 03153 nip = working_space[i + 1][j + 2 * ssizey_ext] / maxch; 03154 nim = working_space[i + 1][j + 1 + 2 * ssizey_ext] / maxch; 03155 for(l = 1; l <= averWindow; l++){ 03156 if(j + l > ymax) 03157 a = working_space[i][ymax + 2 * ssizey_ext] / maxch; 03158 03159 else 03160 a = working_space[i][j + l + 2 * ssizey_ext] / maxch; 03161 b = a - nip; 03162 if(a + nip <= 0) 03163 a = 1; 03164 03165 else 03166 a=TMath::Sqrt(a + nip); 03167 b = b / a; 03168 b = TMath::Exp(b); 03169 spy = spy + b; 03170 if(j - l + 1 < ymin) 03171 a = working_space[i][ymin + 2 * ssizey_ext] / maxch; 03172 03173 else 03174 a = working_space[i][j - l + 1 + 2 * ssizey_ext] / maxch; 03175 b=a-nim; 03176 if(a + nim <= 0) 03177 a = 1; 03178 else 03179 a=TMath::Sqrt(a + nim); 03180 b = b / a; 03181 b = TMath::Exp(b); 03182 smy = smy + b; 03183 } 03184 a = (spx * working_space[i][j + 1] + spy * working_space[i + 1][j]) / (smx + smy); 03185 working_space[i + 1][j + 1] = a; 03186 nom = nom + a; 03187 } 03188 } 03189 for(i = xmin; i <= xmax; i++){ 03190 for(j = ymin; j <= ymax; j++){ 03191 working_space[i][j] = working_space[i][j] / nom; 03192 } 03193 } 03194 for(i = 0; i < ssizex_ext; i++){ 03195 for(j = 0; j < ssizey_ext; j++){ 03196 working_space[i][j + ssizey_ext] = working_space[i][j] * plocha; 03197 working_space[i][2 * ssizey_ext + j] = working_space[i][ssizey_ext + j]; 03198 } 03199 } 03200 } 03201 //deconvolution starts 03202 area = 0; 03203 lhx = -1,lhy = -1; 03204 positx = 0,posity = 0; 03205 maximum = 0; 03206 //generate response matrix 03207 for(i = 0; i < ssizex_ext; i++){ 03208 for(j = 0; j < ssizey_ext; j++){ 03209 lda = (double)i - 3 * sigma; 03210 ldb = (double)j - 3 * sigma; 03211 lda = (lda * lda + ldb * ldb) / (2 * sigma * sigma); 03212 k=(int)(1000 * TMath::Exp(-lda)); 03213 lda = k; 03214 if(lda != 0){ 03215 if((i + 1) > lhx) 03216 lhx = i + 1; 03217 03218 if((j + 1) > lhy) 03219 lhy = j + 1; 03220 } 03221 working_space[i][j] = lda; 03222 area = area + lda; 03223 if(lda > maximum){ 03224 maximum = lda; 03225 positx = i,posity = j; 03226 } 03227 } 03228 } 03229 //read source matrix 03230 for(i = 0;i < ssizex_ext; i++){ 03231 for(j = 0;j < ssizey_ext; j++){ 03232 working_space[i][j + 14 * ssizey_ext] = TMath::Abs(working_space[i][j + ssizey_ext]); 03233 } 03234 } 03235 //calculate matrix b=ht*h 03236 i = lhx - 1; 03237 if(i > ssizex_ext) 03238 i = ssizex_ext; 03239 03240 j = lhy - 1; 03241 if(j>ssizey_ext) 03242 j = ssizey_ext; 03243 03244 i1min = -i,i1max = i; 03245 i2min = -j,i2max = j; 03246 for(i2 = i2min; i2 <= i2max; i2++){ 03247 for(i1 = i1min; i1 <= i1max; i1++){ 03248 ldc = 0; 03249 j2min = -i2; 03250 if(j2min < 0) 03251 j2min = 0; 03252 03253 j2max = lhy - 1 - i2; 03254 if(j2max > lhy - 1) 03255 j2max = lhy - 1; 03256 03257 for(j2 = j2min; j2 <= j2max; j2++){ 03258 j1min = -i1; 03259 if(j1min < 0) 03260 j1min = 0; 03261 03262 j1max = lhx - 1 - i1; 03263 if(j1max > lhx - 1) 03264 j1max = lhx - 1; 03265 03266 for(j1 = j1min; j1 <= j1max; j1++){ 03267 lda = working_space[j1][j2]; 03268 ldb = working_space[i1 + j1][i2 + j2]; 03269 ldc = ldc + lda * ldb; 03270 } 03271 } 03272 k = (i1 + ssizex_ext) / ssizex_ext; 03273 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext] = ldc; 03274 } 03275 } 03276 //calculate at*y and write into p 03277 i = lhx - 1; 03278 if(i > ssizex_ext) 03279 i = ssizex_ext; 03280 03281 j = lhy - 1; 03282 if(j > ssizey_ext) 03283 j = ssizey_ext; 03284 03285 i2min = -j,i2max = ssizey_ext + j - 1; 03286 i1min = -i,i1max = ssizex_ext + i - 1; 03287 for(i2 = i2min; i2 <= i2max; i2++){ 03288 for(i1=i1min;i1<=i1max;i1++){ 03289 ldc=0; 03290 for(j2 = 0; j2 <= (lhy - 1); j2++){ 03291 for(j1 = 0; j1 <= (lhx - 1); j1++){ 03292 k2 = i2 + j2,k1 = i1 + j1; 03293 if(k2 >= 0 && k2 < ssizey_ext && k1 >= 0 && k1 < ssizex_ext){ 03294 lda = working_space[j1][j2]; 03295 ldb = working_space[k1][k2 + 14 * ssizey_ext]; 03296 ldc = ldc + lda * ldb; 03297 } 03298 } 03299 } 03300 k = (i1 + ssizex_ext) / ssizex_ext; 03301 working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext] = ldc; 03302 } 03303 } 03304 //move matrix p 03305 for(i2 = 0; i2 < ssizey_ext; i2++){ 03306 for(i1 = 0; i1 < ssizex_ext; i1++){ 03307 k = (i1 + ssizex_ext) / ssizex_ext; 03308 ldb = working_space[(i1 + ssizex_ext) % ssizex_ext][i2 + ssizey_ext + ssizey_ext + k * 3 * ssizey_ext]; 03309 working_space[i1][i2 + 14 * ssizey_ext] = ldb; 03310 } 03311 } 03312 //initialization in x1 matrix 03313 for(i2 = 0; i2 < ssizey_ext; i2++){ 03314 for(i1 = 0; i1 < ssizex_ext; i1++){ 03315 working_space[i1][i2 + ssizey_ext] = 1; 03316 working_space[i1][i2 + 2 * ssizey_ext] = 0; 03317 } 03318 } 03319 //START OF ITERATIONS 03320 for(lindex = 0; lindex < deconIterations; lindex++){ 03321 for(i2 = 0; i2 < ssizey_ext; i2++){ 03322 for(i1 = 0; i1 < ssizex_ext; i1++){ 03323 lda = working_space[i1][i2 + ssizey_ext]; 03324 ldc = working_space[i1][i2 + 14 * ssizey_ext]; 03325 if(lda > 0.000001 && ldc > 0.000001){ 03326 ldb=0; 03327 j2min=i2; 03328 if(j2min > lhy - 1) 03329 j2min = lhy - 1; 03330 03331 j2min = -j2min; 03332 j2max = ssizey_ext - i2 - 1; 03333 if(j2max > lhy - 1) 03334 j2max = lhy - 1; 03335 03336 j1min = i1; 03337 if(j1min > lhx - 1) 03338 j1min = lhx - 1; 03339 03340 j1min = -j1min; 03341 j1max = ssizex_ext - i1 - 1; 03342 if(j1max > lhx - 1) 03343 j1max = lhx - 1; 03344 03345 for(j2 = j2min; j2 <= j2max; j2++){ 03346 for(j1 = j1min; j1 <= j1max; j1++){ 03347 k = (j1 + ssizex_ext) / ssizex_ext; 03348 ldc = working_space[(j1 + ssizex_ext) % ssizex_ext][j2 + ssizey_ext + 10 * ssizey_ext + k * 2 * ssizey_ext]; 03349 lda = working_space[i1 + j1][i2 + j2 + ssizey_ext]; 03350 ldb = ldb + lda * ldc; 03351 } 03352 } 03353 lda = working_space[i1][i2 + ssizey_ext]; 03354 ldc = working_space[i1][i2 + 14 * ssizey_ext]; 03355 if(ldc * lda != 0 && ldb != 0){ 03356 lda =lda * ldc / ldb; 03357 } 03358 03359 else 03360 lda=0; 03361 working_space[i1][i2 + 2 * ssizey_ext] = lda; 03362 } 03363 } 03364 } 03365 for(i2 = 0; i2 < ssizey_ext; i2++){ 03366 for(i1 = 0; i1 < ssizex_ext; i1++) 03367 working_space[i1][i2 + ssizey_ext] = working_space[i1][i2 + 2 * ssizey_ext]; 03368 } 03369 } 03370 //looking for maximum 03371 maximum=0; 03372 for(i = 0; i < ssizex_ext; i++){ 03373 for(j = 0; j < ssizey_ext; j++){ 03374 working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext] = area * working_space[i][j + ssizey_ext]; 03375 if(maximum < working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext]) 03376 maximum = working_space[(i + positx) % ssizex_ext][(j + posity) % ssizey_ext]; 03377 03378 } 03379 } 03380 //searching for peaks in deconvolved spectrum 03381 for(i = 1; i < ssizex_ext - 1; i++){ 03382 for(j = 1; j < ssizey_ext - 1; j++){ 03383 if(working_space[i][j] > working_space[i - 1][j] && working_space[i][j] > working_space[i - 1][j - 1] && working_space[i][j] > working_space[i][j - 1] && working_space[i][j] > working_space[i + 1][j - 1] && working_space[i][j] > working_space[i + 1][j] && working_space[i][j] > working_space[i + 1][j + 1] && working_space[i][j] > working_space[i][j + 1] && working_space[i][j] > working_space[i - 1][j + 1]){ 03384 if(i >= shift && i < ssizex + shift && j >= shift && j < ssizey + shift){ 03385 if(working_space[i][j] > threshold * maximum / 100.0){ 03386 if(peak_index < fMaxPeaks){ 03387 for(k = i - 1,a = 0,b = 0; k <= i + 1; k++){ 03388 a += (double)(k - shift) * working_space[k][j]; 03389 b += working_space[k][j]; 03390 } 03391 ax=a/b; 03392 if(ax < 0) 03393 ax = 0; 03394 03395 if(ax >= ssizex) 03396 ax = ssizex - 1; 03397 03398 for(k = j - 1,a = 0,b = 0; k <= j + 1; k++){ 03399 a += (double)(k - shift) * working_space[i][k]; 03400 b += working_space[i][k]; 03401 } 03402 ay=a/b; 03403 if(ay < 0) 03404 ay = 0; 03405 03406 if(ay >= ssizey) 03407 ay = ssizey - 1; 03408 03409 if(peak_index == 0){ 03410 fPositionX[0] = ax; 03411 fPositionY[0] = ay; 03412 peak_index = 1; 03413 } 03414 03415 else{ 03416 for(k = 0, priz = 0; k < peak_index && priz == 0; k++){ 03417 if(working_space[shift+(int)(ax+0.5)][15 * ssizey_ext + shift + (int)(ay+0.5)] > working_space[shift+(int)(fPositionX[k]+0.5)][15 * ssizey_ext + shift + (int)(fPositionY[k]+0.5)]) 03418 priz = 1; 03419 } 03420 if(priz == 0){ 03421 if(k < fMaxPeaks){ 03422 fPositionX[k] = ax; 03423 fPositionY[k] = ay; 03424 } 03425 } 03426 03427 else{ 03428 for(l = peak_index; l >= k; l--){ 03429 if(l < fMaxPeaks){ 03430 fPositionX[l] = fPositionX[l - 1]; 03431 fPositionY[l] = fPositionY[l - 1]; 03432 } 03433 } 03434 fPositionX[k - 1] = ax; 03435 fPositionY[k - 1] = ay; 03436 } 03437 if(peak_index < fMaxPeaks) 03438 peak_index += 1; 03439 } 03440 } 03441 } 03442 } 03443 } 03444 } 03445 } 03446 //writing back deconvolved spectrum 03447 for(i = 0; i < ssizex; i++){ 03448 for(j = 0; j < ssizey; j++){ 03449 dest[i][j] = working_space[i + shift][j + shift]; 03450 } 03451 } 03452 i = (int)(4 * sigma + 0.5); 03453 i = 4 * i; 03454 for (j = 0; j < ssizex + i; j++) 03455 delete[]working_space[j]; 03456 delete[]working_space; 03457 fNPeaks = peak_index; 03458 return fNPeaks; 03459 } 03460 03461 03462 // STATIC functions (called by TH1) 03463 03464 //_______________________________________________________________________________ 03465 Int_t TSpectrum2::StaticSearch(const TH1 *hist, Double_t sigma, Option_t *option, Double_t threshold) 03466 { 03467 //static function, interface to TSpectrum2::Search 03468 03469 TSpectrum2 s; 03470 return s.Search(hist,sigma,option,threshold); 03471 } 03472 03473 //_______________________________________________________________________________ 03474 TH1 *TSpectrum2::StaticBackground(const TH1 *hist,Int_t niter, Option_t *option) 03475 { 03476 //static function, interface to TSpectrum2::Background 03477 03478 TSpectrum2 s; 03479 return s.Background(hist,niter,option); 03480 }
1.5.1