Dear all,
TSpectrum *s = new TSpectrum(npeaks);
Int_t nfound = s->Search(peakno,0.0,"",0.1);
Int_t peaks;
auto *xpeaks = s->GetPositionX();
for (Int_t p=0;p<nfound;p++ )
{
Float_t xp =xpeaks[p];
std::cout << "xp " << xp << std::endl;
Int_t bin =peakno->GetXaxis()->FindBin(xp);
std::cout << "bin = " << bin << std::endl;
Float_t yp=peakno->GetBinContent(bin);
std::cout << "yp = " <<yp << std::endl;
std::cout << " p = " << p << "nfound = " << nfound << std::endl;
xpeaks++;
}
The output is like
yp = 1
p = 4nfound = 10
xp 0
bin = 1
yp = 0
p = 5nfound = 10
xp 2
bin = 3
yp = 0
p = 6nfound = 10
xp 1
bin = 2
yp = 0
p = 7nfound = 10
yp = 0
p = 3nfound = 4
go on 1000
xp 6789.5
bin = 6790
yp = 39
p = 0nfound = 3
xp 6814.5
bin = 6815
yp = 6
p = 1nfound = 3
xp 0
bin = 1
yp = 0
p = 2nfound = 3
yp = 0
p = 1nfound = 2
go on 1000
xp 6789.5
bin = 6790
yp = 119
p = 0nfound = 2
xp 0
bin = 1
yp = 0
p = 1nfound = 2
here the value of nfound is changing. I need to print out the final value that is nfound =2 at the end of the loop. The code is printing the value of nfound throughout the loop but i need only the final value at the end of the loop. How do I can do that ? Could anybody help ?
Thank you
yus
November 21, 2023, 12:12pm
2
Hi,
nfound does not change inside the loop, so the only reason you see it changing is because you apparently have several spectra for which you are trying to search the peaks.
Also, could you please share the exact code you are using? Why do you have the linesgo on 1000
Here is the code ,
#include <iostream>
#include <string>
#include <sstream>
#include <fstream>
#include <unistd.h>
#include <stdlib.h>
#include <algorithm>
#include "struct.h"
#include "ReadGDML.h"
#include "ReadEDepSim.h"
#include "ReadOptical.h"
#include "ReadDigit.h"
#include "CameraImage.h"
#include "FastFitter.h"
#include "LensDisplayer.h"
#include "LineIntersector.h"
#include "StructDefinition.h"
#include "TrackMatcher.h"
#include "Volume3DReconstructer.h"
#include "EnergyEstimator.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TApplication.h"
#include "TF1.h"
#include "TH2I.h"
#include "TCanvas.h"
#include "TVector3.h"
#include "TH1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TChain.h"
#include "TVirtualFitter.h"
#include "Riostream.h"
Int_t npeaks=10 ;
Int_t final_nfound;
Int_t final_nfound1;
Int_t A;
Int_t B;
Int_t t;
float taufast=6e-9;
float tauslow=1.6e-6;
void usage() {
std::cout << "Input parameters needed: " << std::endl;
std::cout << "-i 'root file with spill data '" << std::endl;
std::cout << std::endl;
}
int main(int argc, char **argv) {
std::string inFilename="";
std::string outFilename="";
int iflag=0, oflag=0; //mandatory
int opt;
outFilename = "interaction.root";
while( (opt = getopt(argc, argv, "i:o")) != -1) {
switch(opt) {
case 'o':
if(optarg) outFilename = optarg;
oflag = 1;
break;
case 'i':
if(optarg) inFilename = optarg;
iflag = 1;
break;
default:
usage();
return 1;
}
}
if(!iflag) {
std::cout << "ERROR: missing REQUIRED parameters!" << std::endl;
usage();
return 1;
}
std::cout << "Sensor file: " << inFilename << std::endl;
std::cout << "Output file: " << outFilename << std::endl;
std::cout << "*****************************************************" << std::endl;
std::cout << "Reading data" << std::endl;
// open file
TFile *m_f = new TFile(inFilename.c_str(), "READ");
if (!m_f) {
std::cout << "ERROR: missing the file!" << std::endl;
}
Int_t k;
TChain *tree = new TChain();
for (int k=1;k<39;k++){
tree->Add(Form("./%s/CAM_%d",inFilename.c_str(),k));
}
// TTree *tree= (TTree *)m_f->Get("CAM_4");
if (!tree) {
std::cout << "ERROR: missing the event trees!" << std::endl;
}
Int_t idEvent;
std::vector<double> *energy=0 ;
std::vector<double> *time=0 ;
tree->SetBranchAddress("time", &time);
tree->SetBranchAddress("idEvent",&idEvent);
tree->SetBranchAddress("energy",&energy);
TH1F *peakno = new TH1F("peakno","peakno",20000,0,20000);
TH1F *intno = new TH1F("intno","intno",100,0,100);
int N = tree->GetEntries();
// intno->Fill(final_nfound);
std::cout << "N= " << N << std::endl;
// std::cout << "go on 50 " << std::endl;
for (int i=0; i<N; i++)
{
tree->GetEntry(i);
ULong_t nsig = time->size();
double *sig = time->data();
// std::cout << " go on 70 " << std::endl;
double *esig = energy->data();
if (idEvent==0){
for (ULong_t j=0; j<nsig; j++){
peakno->Fill(sig[j]);
}
// std::cout << " go on 500 " << std::endl;
TSpectrum *s = new TSpectrum(npeaks);
Int_t nfound = s->Search(peakno,0.0,"",0.1);
Int_t peaks;
auto *xpeaks = s->GetPositionX();
// std::cout << " go on 700 " << std::endl;
for (Int_t p=0;p<nfound;p++ )
{
Float_t xp =xpeaks[p];
// std::cout << " go on 150" << std::endl;
std::cout << "xp " << xp << std::endl;
Int_t bin =peakno->GetXaxis()->FindBin(xp);
std::cout << "bin = " << bin << std::endl;
Float_t yp=peakno->GetBinContent(bin);
std::cout << "yp = " <<yp << std::endl;
std::cout << " p = " << p << "nfound = " << nfound << std::endl;
xpeaks++;
// intno->Fill(nfound);
}
// intno->Fill(nfound);
std::cout << "final_nfound=" << final_nfound1 << std::endl;
std::cout << " go on 1000 " << std::endl;
}
// std::cout << " go on 1500 " << std::endl;
}
// std::cout << "go on 2000 " << std::endl;
TFile *outputFile = new TFile(outFilename.c_str(),"RECREATE");
std::cout << "Output file written to " << outFilename << std::endl;
peakno->Write();
intno->Write();
// std::cout << " go on 2500 " << std::endl;
}
go on 1000 these are just to check the sequence of the code.
thanks
yus
November 21, 2023, 12:44pm
4
All right, thanks. Can it so happen that you have multiple entries in your TTree with idEvent==0? Or is it supposed to be unique?
Yes, idEvent == 0 supposed to be unique .
yus
November 21, 2023, 12:56pm
6
You use a TChain, meaning that you use multiple input files. Can it so happen that idEvent==0 is true for more than one entry in all these input files combined?
Yes, there are total 38 trees . All trees contain a spectrum . I added all the trees to get a single spectra. Then using peaksearch trying to find the number of peaks from the total spectra.
idevent==0 is true for every trees.
yus
November 21, 2023, 1:34pm
8
Which means you may have an entry with idEvent==0 in the first file, then in the second file, and so on? And the spectra will always be different and have different numbers of found peaks, no?
yes. But the final result gives the peak no of the added spectra.
I think the spectrums are not adding correctly … that is the problem.
yus
November 21, 2023, 1:44pm
11
I don’t see where this is done. You have a loop over entries in all your files:
int N = tree->GetEntries();
...
for (int i=0; i<N; i++)
{
tree->GetEntry(i);
...
}
Where do you add the spectra across all your files?
When I am adding using TChain all the trees , all the spectrums will be added automatically. Right ?
In the final spectrum the counts are also increased. The problem is while filling the histogram with nfound it is filling each value of nfound . Not with the final value. Suppose finally after peak searching i get 2 no of peaks , i want to fill the histogram with this no 2 with count 1 .
in this final spectra the two peaks are clearly identified with good counts.
yus
November 21, 2023, 2:03pm
15
I don’t think so. By using TChain you add TTree entries from multiple files, that’s it.
Then How do i can do that ? I can open only one tree and read the spectra . But i need to add all the trees to get the added spectra. i can not open all the trees together.
couet
November 21, 2023, 2:13pm
17
I have indented you code for a better readability:
#include <iostream>
#include <string>
#include <sstream>
#include <fstream>
#include <unistd.h>
#include <stdlib.h>
#include <algorithm>
#include "struct.h"
#include "ReadGDML.h"
#include "ReadEDepSim.h"
#include "ReadOptical.h"
#include "ReadDigit.h"
#include "CameraImage.h"
#include "FastFitter.h"
#include "LensDisplayer.h"
#include "LineIntersector.h"
#include "StructDefinition.h"
#include "TrackMatcher.h"
#include "Volume3DReconstructer.h"
#include "EnergyEstimator.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TApplication.h"
#include "TF1.h"
#include "TH2I.h"
#include "TCanvas.h"
#include "TVector3.h"
#include "TH1.h"
#include "TRandom.h"
#include "TSpectrum.h"
#include "TChain.h"
#include "TVirtualFitter.h"
#include "Riostream.h"
Int_t npeaks=10 ;
Int_t final_nfound;
Int_t final_nfound1;
Int_t A;
Int_t B;
Int_t t;
float taufast=6e-9;
float tauslow=1.6e-6;
void usage() {
std::cout << "Input parameters needed: " << std::endl;
std::cout << "-i 'root file with spill data '" << std::endl;
std::cout << std::endl;
}
int main(int argc, char **argv) {
std::string inFilename="";
std::string outFilename="";
int iflag=0, oflag=0; //mandatory
int opt;
outFilename = "interaction.root";
while( (opt = getopt(argc, argv, "i:o")) != -1) {
switch(opt) {
case 'o':
if(optarg) outFilename = optarg;
oflag = 1;
break;
case 'i':
if(optarg) inFilename = optarg;
iflag = 1;
break;
default:
usage();
return 1;
}
}
if (!iflag) {
std::cout << "ERROR: missing REQUIRED parameters!" << std::endl;
usage();
return 1;
}
std::cout << "Sensor file: " << inFilename << std::endl;
std::cout << "Output file: " << outFilename << std::endl;
std::cout << "*****************************************************" << std::endl;
std::cout << "Reading data" << std::endl;
// open file
TFile *m_f = new TFile(inFilename.c_str(), "READ");
if (!m_f) {
std::cout << "ERROR: missing the file!" << std::endl;
return 1;
}
Int_t k;
TChain *tree = new TChain();
for (int k=1;k<39;k++) {
tree->Add(Form("./%s/CAM_%d",inFilename.c_str(),k));
}
if (!tree) {
std::cout << "ERROR: missing the event trees!" << std::endl;
}
Int_t idEvent;
std::vector<double> *energy=0 ;
std::vector<double> *time=0 ;
tree->SetBranchAddress("time", &time);
tree->SetBranchAddress("idEvent",&idEvent);
tree->SetBranchAddress("energy",&energy);
TH1F *peakno = new TH1F("peakno","peakno",20000,0,20000);
TH1F *intno = new TH1F("intno","intno",100,0,100);
int N = tree->GetEntries();
std::cout << "N= " << N << std::endl;
for (int i=0; i<N; i++) {
tree->GetEntry(i);
ULong_t nsig = time->size();
double *sig = time->data();
double *esig = energy->data();
if (idEvent==0){
for (ULong_t j=0; j<nsig; j++) peakno->Fill(sig[j]);
TSpectrum *s = new TSpectrum(npeaks);
Int_t nfound = s->Search(peakno,0.0,"",0.1);
Int_t peaks;
auto *xpeaks = s->GetPositionX();
// std::cout << " go on 700 " << std::endl;
for (Int_t p=0;p<nfound;p++ ) {
Float_t xp =xpeaks[p];
std::cout << "xp " << xp << std::endl;
Int_t bin =peakno->GetXaxis()->FindBin(xp);
std::cout << "bin = " << bin << std::endl;
Float_t yp=peakno->GetBinContent(bin);
std::cout << "yp = " <<yp << std::endl;
std::cout << " p = " << p << "nfound = " << nfound << std::endl;
xpeaks++;
}
std::cout << "final_nfound=" << final_nfound1 << std::endl;
}
}
TFile *outputFile = new TFile(outFilename.c_str(),"RECREATE");
std::cout << "Output file written to " << outFilename << std::endl;
peakno->Write();
intno->Write();
}
There is something I do not undertstand. You print the the final value of nfound with the line:
std::cout << "final_nfound=" << final_nfound1 << std::endl;
but you never set the variable final_nfound1
I am sorry . This line was a mistake . I was checking something and forget to comment out.
Is it possible to open more than 1 trees containing spectra and make a addition of them to get a single spectra ?
couet
November 21, 2023, 2:29pm
20
If you start from existing histograms you can add histograms (see doc).
If you start from the trees and chain them, then the histogram you will create will embrase the whole data set and the peaks contained in the chained-trees will be reflected in that histogram.
