【ROOT from CERN】——TSpectrum2类与二维寻峰

已于 2022-03-05 17:42:58 修改
阅读量1.6k 收藏 10
点赞数 3
分类专栏： ROOT 文章标签： c++
于 2021-11-18 00:39:53 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.youkuaiyun.com/qq_20604319/article/details/121388282
版权

TSpectrum2类是二维谱类，其中的方法可以寻找在二维直方图（即散点图）内的峰个数。本篇文章会介绍部分该类的基本用法和原理。

一、基本简介

根据官网介绍，TSpectrum类以及TSpectrum2类是遗留代码，其在以后的ROOT版本中已经不再继续维护，但仍可使用。该类在《ROOTUsersGuide》也未列出，因此只能查看官网的电子文档来获取相关信息。以下引自官网声明：

Legacy Code

The Spectrum package is a legacy interface: it is not recommended to use it in new code. There will be no bug fixes nor new developments.

关于各种电子文档的查阅，各位读者可以尝试在ROOT Manual - ROOT该链接中查找，还请读者自己探索。而本文提及的TSpectrum类和TSpectrum2类可以在ROOT: Advanced spectra processing classes.之中找到。网站中使用的范例代码和数据，各位均可以在$ROOTSYS/tutorials/下找到。

二、寻峰代码

1、方法SearchHighRes的信息

Int_t TSpectrum2::SearchHighRes 	( 	Double_t **  	source,
		Double_t **  	dest,
		Int_t  	ssizex,
		Int_t  	ssizey,
		Double_t  	sigma,
		Double_t  	threshold,
		Bool_t  	backgroundRemove,
		Int_t  	deconIterations,
		Bool_t  	markov,
		Int_t  	averWindow 
	)

这是二维寻峰的核心方法的函数原型，简单介绍一下各个参数：（翻译可能不准确，有些参数还未完全弄懂）

  • source - 指向源谱矩阵的指针
  • dest - 指向解反卷积谱矩阵的指针
  • ssizex - 源谱x轴长度
  • ssizey - 源谱y轴长度
  • sigma - 搜索峰值的sigma，详细信息请参考手册（关于手册中的补充：峰值中心的值 [ i ] 减去两个对称的channel长度（channels i-3*sigma, i+3*sigma）的平均值必须大于阈值threshold。 否则峰值将被忽略。可以近似理解为正态分布中的sigma。）
  • threshold - 所选峰值的阈值，幅度小于threshold*highest_peak/100的峰将被忽略，参见手册
  • backgroundRemove - 逻辑变量，设置是否需要在反卷积之前去除背景
  • deconIterations - 反卷积操作中的迭代次数
  • markov - 逻辑变量，如果为真，首先将源谱替换为使用马尔可夫链方法计算的新谱。
  • averWindow - 搜索峰值的平均窗口，详见手册（仅适用于马尔可夫链方法）

为了防止歧义，把原英文描述贴在下面：

  • source-pointer to the matrix of source spectrum
  • dest-pointer to the matrix of resulting deconvolved spectrum
  • ssizex-x length of source spectrum
  • ssizey-y length of source spectrum
  • sigma-sigma of searched peaks, for details we refer to manual
  • threshold-threshold value in % for selected peaks, peaks with amplitude less than threshold*highest_peak/100 are ignored, see manual
  • backgroundRemove-logical variable, set if the removal of background before deconvolution is desired
  • deconIterations-number of iterations in deconvolution operation
  • markov-logical variable, if it is true, first the source spectrum is replaced by new spectrum calculated using Markov chains method.
  • averWindow-averaging window of searched peaks, for details we refer to manual (applies only for Markov method)

2、分类介绍 

以下介绍该代码使用示例为上述网址内的，Example 8 - Src.C。我们分三个部分介绍。

#include <TSpectrum2.h>
 
void Src() {
   Int_t i, j, nfound;
   const Int_t nbinsx = 64;
   const Int_t nbinsy = 64;
   Double_t xmin = 0;
   Double_t xmax = (Double_t)nbinsx;
   Double_t ymin = 0;
   Double_t ymax = (Double_t)nbinsy;
   Double_t** source = new Double_t*[nbinsx];
   for (i=0;i<nbinsx;i++)
      source[i]=new Double_t[nbinsy];
   Double_t** dest = new Double_t*[nbinsx];
   for (i=0;i<nbinsx;i++)
      dest[i]=new Double_t[nbinsy];
   TString dir  = gROOT->GetTutorialDir();
   TString file = dir+"/spectrum/TSpectrum2.root";
   TFile *f     = new TFile(file.Data());
   gStyle->SetOptStat(0);
   auto search = (TH2F*) f->Get("search4");
   auto *s = new TSpectrum2();
   for (i = 0; i < nbinsx; i++){
      for (j = 0; j < nbinsy; j++){
         source[i][j] = search->GetBinContent(i + 1,j + 1);
      }
   }

/**************************************************************************************/

   nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);

/**************************************************************************************/

   printf("Found %d candidate peaks\n",nfound);
   Double_t *PositionX = s->GetPositionX();
   Double_t *PositionY = s->GetPositionY();
   search->Draw("COL");
   auto m = new TMarker();
   m->SetMarkerStyle(23);
   m->SetMarkerColor(kRed);
   for (i=0;i<nfound;i++) {
      printf("posx= %d, posy= %d, value=%d\n",(Int_t)(PositionX[i]+0.5), (Int_t)(PositionY[i]+0.5),
      (Int_t)source[(Int_t)(PositionX[i]+0.5)][(Int_t)(PositionY[i]+0.5)]);
      m->DrawMarker(PositionX[i],PositionY[i]);
   }
}

1、声明部分

Int_t i, j, nfound;
   const Int_t nbinsx = 64;//it is the same as source spectrum
   const Int_t nbinsy = 64;//it is the same as source spectrum
   Double_t xmin = 0;
   Double_t xmax = (Double_t)nbinsx;
   Double_t ymin = 0;
   Double_t ymax = (Double_t)nbinsy;
   Double_t** source = new Double_t*[nbinsx];
   for (i=0;i<nbinsx;i++)
      source[i]=new Double_t[nbinsy];
   Double_t** dest = new Double_t*[nbinsx];
   for (i=0;i<nbinsx;i++)
      dest[i]=new Double_t[nbinsy];
   TString dir  = gROOT->GetTutorialDir();
   TString file = dir+"/spectrum/TSpectrum2.root";
   TFile *f     = new TFile(file.Data());//you can change this root filename
   gStyle->SetOptStat(0);
   auto search = (TH2F*) f->Get("search4");//you can change this histogram name
   auto *s = new TSpectrum2();
   for (i = 0; i < nbinsx; i++){
      for (j = 0; j < nbinsy; j++){
         source[i][j] = search->GetBinContent(i + 1,j + 1);
      }
   }

SearchHighRes()中的source, dest, nbinsx, nbinsy在此处被定义，其中大部分无需我们修改。我们可以操作的参数只有nbinsx和nbinsy，其与源谱应保持一致。它的含义是，对粒子谱的x轴的0-64和y轴的0-64的区域内进行寻峰。注意：1、nbinsx和nbinsy的上限是100。也就是说最大搜索区间是0-100。2、强调其搜索区间是大于零的，只能在第一象限。（这是笔者多次尝试后得到的结果，如果读者发现错误，还请指正。）3、该方法source指向一个二维直方图TH2F类（对于该示例指向search4）。故如果源谱的坐标轴范围超出了100，则需要对填充二维直方图的各个粒子坐标进行数学变换，使其压缩到0-100的范围内。因此源谱的轴范围上限是多少，两个参数的值即为多少。

除注释外，以上结构无需变动。读者可以到前文提到的网址或root的文件夹下查看其他示例来总结不变的代码框架是什么，实际上就是共有的代码部分。

2、核心部分

nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);

该核心部分的各参数含义已经在上文提及过了。出去前面声明部分已经定义了一部分参数，现在在该语句中我们可以操作的参数有sigma，threshold，backgroundRemove，deconlterations，markov，averWindow。

其中backgroundRemove与deconlterations共同使用。前者为KTRUE时，后者的数值才有意义。该迭代次数越多，可以认为对背景的去除越强，同时可能让一些微小的峰变得显著。

其中markov与averWindow共同使用。前者为KTRUE时，后者的数值才有意义。开启马尔科夫链算法后，源谱会被一个更加平滑的谱所替代，可以有效降低源谱的躁动和“参差程度”。平均窗口数值越大，生成的新谱越光滑。

如上图所示，红色线为马尔科夫链算法生成的新谱（此示例为一维谱），其波动远小于绿色的源谱，变得更加平滑。图片引自官方出具的手册《Spectrum》。

该函数返回寻峰的个数（Int_t类型），此代码用指针s访问其返回值并赋值给nfound。对于该函数的一组参数，其结果可能只对具有某类特征的峰显著性和正确性，因此在面对不同情况的源谱时，需要先调参以确保结果正确。

3、绘图部分

   printf("Found %d candidate peaks\n",nfound);
   Double_t *PositionX = s->GetPositionX();
   Double_t *PositionY = s->GetPositionY();
   search->Draw("COL");
   auto m = new TMarker();
   m->SetMarkerStyle(23);
   m->SetMarkerColor(kRed);
   for (i=0;i<nfound;i++) {
      printf("posx= %d, posy= %d, value=%d\n",(Int_t)(PositionX[i]+0.5), (Int_t)(PositionY[i]+0.5),
      (Int_t)source[(Int_t)(PositionX[i]+0.5)][(Int_t)(PositionY[i]+0.5)]);
      m->DrawMarker(PositionX[i],PositionY[i]);

设置了一些绘图的输出设置。无需改动。如果只想知道SearchHighRes()的返回值，则可以把该模块全部删除，仅输出nfound即可。可以减小图形开销。当然，更加直观的图片结果更有利于程序员来判断代码的正误，进行调试。因此在对SearchHighRes()调参时，最好保留该模块。

以下两张图分别是源谱和寻峰后的生成谱，其中红色三角标识是寻到的峰位。

 

有关于TSpectrum2类的寻峰代码介绍到这。简单介绍，如需学习，还请各位查阅前文提到的官方网站和官方出具的手册进行学习

【资料】

1、ROOT官网——ROOT: analyzing petabytes of data, scientifically. - ROOT

2、ROOT官方扩展手册——《Spectrum》

如有错误请指正。

java二维数组算法_关于Java：在二维数组中查找峰的算法
weixin_33426157的博客
02-28 546
假设我在java int[][] array中有一个2D累加器数组。该数组如下所示：(x和z轴表示数组中的索引，y轴表示值-这些是int[56][56]的图像，其值介于0?4500之间)要么我需要做的是在阵列中找到峰-第一个中有2个峰，第二个阵列中有8个峰。这些峰值始终是"明显的"(峰值之间始终存在间隙)，但是它们不必像这些图像一样相似，它们或多或少是随机的-这些图像不是基于真实数据，而只是基于样...
Matlab 求二维矩阵的2个或者更多个峰值
最新发布
龟速派的博客
10-12 3259
求一个波峰的方式我就不多说了,find函数就可解决，求两个三个或者更多个该怎么解决？
Matlab二维矩阵(三维图像)寻找波峰
03-24
用于matlab寻找三维图像或二维矩阵的波峰，详细使用方法可在 https://github.com/JamesHoi/Matlab/tree/master/IMMC/Findpeaks 查看。此资源免费下载，我的更多资源可在 https://github.com/JamesHoi?tab=repositories 下载
山脉数组的峰顶索引(二分查找法)
weixin_42738495的博客
01-15 944
我们把符合下列属性的数组 A 称作山脉： A.length >= 3 存在 0 < i < A.length - 1 使得A[0] < A[1] < … A[i-1] < A[i] > A[i+1] > … > A[A.length - 1] 给定一个确定为山脉的数组，返回任何满足 A[0] < A[1] < … A[i-1] &l...
谱峰搜索传播算子（python）
wudongming2017的博客
12-23 1439
使用谱峰搜索传播算子的方法进行DOA估计
二维搜索峰值matlab代码-real_time_vibrotactile_feedback_3D_trajectories:请参阅下面此页面
06-17
二维搜索最大matlab代码3D 轨迹的实时振动触觉反馈 通过将实时值分别与目标路径和目标速度曲线进行比较，独立地为两个运动分量（位置和速度）提供振动触觉反馈。 复合路径是通过沿轨迹引入通路点来创建的，目标路径是为每个人动态构建的，以考虑诸如高度、臂长等变量。目标速度剖面是使用最小加加速度模型构建的。 1速度被用作量化运动平滑度的代理，它代表运动最优性的度量。 该代码还包含用于 3D 轨迹的运动分解算法。 2任何轨迹都可以看作是最小加加速度子运动的线性组合。 通过最小化依赖于参数的重构切向速度分布与运动的实际切向速度分布之间的误差来确定子运动。 先决条件 Polhemus Liberty 运动追踪器 MATLAB R2007b 触觉执行器 收集数据 运行算法的步骤详述如下： 为实验和对照条件两者： 通过指定沿轨迹应出现过Kong点的分数点（m1、m2 等），在starting_point.m 中设置所需的过Kong点数。 以弧度（theta1、theta2 等）为单位设置角度，在这些角度应在starting_point.m （绕z 轴旋转）中引入通Kong点。 运行fixed_coll
CERN ROOT 安装指导
01-09
 CERN的root的是一个很棒（具体有多棒我也不知道）的数据处理软件，经常需要用到。前两次的安装分别是在两位不同的师兄指导下完成的，过程记得很模糊，且linux基础很差，很多指令不甚理解，可以说是一头雾水。  第...
root_music.rar_ROOT_root MUSIC_root_music_root——music
09-14
"ROOT"在这里通常指的是CERN开发的一个数据分析和可视化软件，广泛应用于高能物理实验的数据处理。"MUSIC"可能是这个程序的名字，或者与音乐数据的处理方式有关。 【描述】"一段关于root_music的仿真程序，程序中对...
CERN ROOT.pptx
07-23
以上文档内容包含链接和中间个人安装中出现的问题和解决方法与测试，且附带一些常见小问题的解决方法。测试安装环境为centos7,
网址：root.cern
02-10
来自root-project / web存储库本身的分支的拉取请求部署在https://root.cern/<PRNUMBER> 。 来自此存储库的派生的拉取请求将不会部署。 请注意，只要部署了新版本的网站，就会删除所有PR部署。 只需重新运行PR的...
7.谱峰搜索传播算子DOA估计算法MATLAB程序.rar_DOA_DOA估计_传播算子_搜索谱峰_谱峰搜索
07-15
对谱峰搜索传播算子DOA估计算法MATLAB程序
DOA估计的谱峰搜索
04-06
关于给定信源数的谱峰搜索，适合精确获取谱峰值
谱峰搜索传播算子算法(PM)实现波达角估计
11-26
使用谱峰搜索传播算子算法(PM)实现波达角估计，基于matlab环境
matlab寻峰代码-MATLAB-Signal-Analysis-Fourier-Transform:用于查找未知信号频率的MATLABFF
06-04
matlab寻峰代码MATLAB-信号分析-傅立叶变换 让我们从傅立叶先生的诺贝尔名言开始：“每个信号都可以表示为加权复指数的总和。” MATLAB FFT() 函数用于寻找未知信号频率 我们已经取了一个随机信号（未知频率），用于 1500 Hz 频率的仿真。 它以比信号频率大 10 倍的频率进行采样。 代码运行并将未知信号数据点作为输入（您必须知道在 matlab 中一切都是数字/离散的，您不能在 matlab 中创建真实世界的信号/模拟信号）。 并在使用 FFT 命令后，将数据点转换为复数值。 实部代表该样本的幅度，虚部代表该样本的相位。 这些样本只不过是不同的频率。 我们如何解释这些频率？ 现在进入信号的重要参数RESOLUTION，它以采样频率/采样数给出。 FFT 中的每个采样点代表一个分辨率(R)。 例如，第一个样本位于第 2 个位置（结果图像中的第二个数字），因为我们在 0 频率处没有任何东西（虚概念）。 信号分析：要找出哪个频率分量对给定信号的贡献最大，请找到第 16 个位置的峰值，因此未知信号的最大频率为 (15 * RESOLUTION) 在这里，我们得到 ~100
matlab寻峰代码-scChIC-seq:scChIC-seq
06-04
matlab寻峰代码单细胞染色质免疫切割测序 (scChIC-Seq) 以分析组蛋白修饰 信息 该存储库包括许多部分，描述了论文中的分析。 请注意，图编号可能与论文编号不同。 这里的分析顺序与论文中的顺序不同。 欢迎任何意见和问题。 A.下载文件 从github下载文件并将目录更改为CRK [wlku@matrix] git clone https://github.com/wailimku/scChIC-seq.git [wlku@matrix] mv scChIC-seq CRK [wlku@matrix] cd CRK [wlku@matrix CRK] ls [wlku@matrix CRK] AdvancedColormap.m data Figures GSE105012 README.md src violin.m violinplot.m 从 GEO 网站下载 GSE105012_RAW.tar 并将其保存到文件夹GSE105012 [wlku@matrix CRK] cd GSE105012 [wlku@matrix GSE105012] ls [wlku@matri
ROOT.jl：Julia的CERN ROOT兼容性
02-03
**ROOT.jl：Julia与CERN ROOT的互操作性** **标题解析：** "ROOT.jl：Julia的CERN ROOT兼容性" 指的是一个名为 `ROOT.jl` 的Julia语言包，其目标是实现Julia与CERN（欧洲核子研究组织）开发的数据分析和科学计算...
[MIT6.006 算法导论] 1. Peak Finding 寻峰
weixin_44179561的博客
04-11 3611
一、1D情况 假设有一个如下图的一维数组，格子下的数字代表它们的索引位置，格子内的字母代表该位置内的数值。 峰值的定义：当且仅当 b≥a 且 b≥c 时，位置2为峰值。（如果位置位于数组两个边界，则只需要比较另一个方向。如当 a≥b 时，位置1为峰值。） 寻峰：如果峰值存在，找到峰值所在的位置。（注意，由于定义时使用了大于等于，在此定义下数组内必定至少存在一个位置满足峰值的定义，因此此时不需要考虑峰值不存在的情况） 1、遍历 最简单且最直接的峰值寻找方式就是遍历，遍历数组内的每一个位置，若该位置内的元素满
在matlab中如何进行谱峰搜索,谱峰搜索的神经网络实现
weixin_42112658的博客
03-23 515
[1]Schmidt R O. IEEE Trans. on AP, 1986, 34(3):276-280.[2]Kumaresan R, Tufts D W. IEEE Trans. on AES, 1983, AES-19(1):134-139.[3]Stoica P, Nehorai A. IEEE Trans. on ASSP, 1989, ASSP-37(5):720-741.[4]S...
二维峰值查找
diker 的博客
03-18 4769
问题描述： 给定一个二维数组，求出其中任意一个峰值。峰值定义为比上下左右都大或相等的元素值（规定数组外无穷小），即 arr[i][j]arr[i][j]arr[i][j] 为峰值当 arr[i][j]arr[i][j]arr[i][j] ≥\geq≥ arr[i−1][j]arr[i-1][j]arr[i−1][j] , arr[i][j]arr[i][j]arr[i][j] ≥\geq≥ arr...
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值