ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided
Hello
I am trying to locate the right most peak from the histogram. However I am having problem where the parameters I’ve found from using GetPositionX() are printing out gibberish. The problem starts from line 27, where I set the value of xpeaks.
struct hit_eff_vs_energy
{
float energy;
float sigma;
float quartz_hits;
float lyso_hits;
float hit_efficiency;
float total_PE;
float PE_eff;
float PE_eff_err;
};
hit_eff_vs_energy eff;
//photon beam energy from run 10942 to 10983
const float Conv[50] = {21.06, 22.11, 23.80, 25.85, 15.43, //10942-10946
17.52, 18.21, 20.31, 26.58, 28.60, //10947-10951
29.30, 31.30, 32.00, 34.00, 34.60, //10952-10956
39.20, 41.80, 44.30, 46.80, 49.30, //10957-10961
51.70, 54.10, 56.4};//10962-10964
//10965 pedestal run:
//Quartz: 166
//TR: 228
//TL: 136
//LYSO: 146
//local calibration: 2.29pC/PE where 0.2pC/channel
float quartz_pedestal = 165.2;
float OnePE = 11.45;
float quartz_ped_cut = quartz_pedestal+OnePE/2; //pedestal =165.2, signle PE = 165.2+10.4 = 175.6, cut at half PE 170
int sigma_cut = 3; // number of sigma used in cuuting on LYSO photopeak
void AnalyzeSingle2(int run)
{
// if (run>=10968)
// pedestal=167;
// TCanvas *c1 = new TCanvas("c1","Mainz Photon Beam Data",800,900);
TCanvas *c1 = new TCanvas();
c1->Divide(2,2);
TCanvas *c2 = new TCanvas("c2","Mainz Photon Beam Data",600,480);
TFile *f_data = TFile::Open(Form("run_%d.root",run));
TTree *dTree = (TTree*)f_data->Get("dataTree");
TH1D *h_lyso = new TH1D("h_lyso","",4096,0,4096);
TSpectrum *s = new TSpectrum(3);
TH1D *h4 = new TH1D("h4","",4096,0,4096);
TH1D *h3 = new TH1D("h3","",4096,0,4096);
c1->cd(1);
TCanvas *c11= new TCanvas();
dTree->Draw("LYSOChn>>h_lyso","");
h_lyso->SetLineColor(4);
//Use TSpectrum to find the peak candidates
float temp_sigma= 20;
// if(Conv[run-10942]>55)
// temp_sigma= 100;
//options: Search(histogram_pointer,sigma of searched peaks,options,threshold)
Int_t nfound = s->Search(h_lyso,temp_sigma,"new",0.03);
printf("Found %d candidate peaks to fitn\n",nfound);
// TH1 *hb= s->Background(h_lyso);
// TH1 *hb= s->Background(h_lyso,5,"BackOrder4");
// hb->Draw("same");
// c11->SetLogy();
//c11->SaveAs("testgraph1.jpg");
//Locate the right most peak
Double_t *xpeaks = s->GetPositionX();
cout<<"Position x: "<<s->GetPositionX()<<endl;
cout<<"xpeaks: "<<xpeaks[1]<<endl;
Double_t *ypeaks = s->GetPositionY();
Double_t Xtemp = 0;
cout<<"xpeaks: "<<sizeof(*xpeaks)<<endl;
for (int p=0;p<(nfound);p++)
{
Double_t xp = xpeaks[p];
if(Xtemp<xp)
Xtemp = xp;
Int_t bin = h_lyso->GetXaxis()->FindBin(xp);
Float_t yp = h_lyso->GetBinContent(bin);
cout<<"xp: "<<xp<<endl;
cout<<"Xtemp: "<<Xtemp<<endl;
cout<<"bin: "<<bin<<endl;
cout<<"yp: "<<yp<<endl;
cout<<"\n"<<endl;
}
printf("Now fitting main photopeak.");
Double_t peak_mean = Xtemp;
cout<<"peak_mean: "<<peak_mean<<endl;
// cout<<"peak_mean: "<<peak_mean<<endl;
Double_t peak_sigma = 2;
TF1 *fitF = new TF1("fitF","gaus(0)",peak_mean-0.3*peak_sigma,peak_mean+1.2*peak_sigma);
for(int l=0; l<5; l++)
{
h_lyso->Fit(fitF,"","",peak_mean-0.3*peak_sigma,peak_mean+1.2*peak_sigma);
peak_mean = fitF->GetParameter(1);
peak_sigma = fitF->GetParameter(2);
}
if (fitF)= c11->Update();
c11->SaveAs("testgraph11.jpg");
eff.energy = peak_mean - 125; // pedestal center of LYSO channel
eff.sigma = peak_sigma;
cout<<"Mean: "<<peak_mean<<", Sigma: "<<peak_sigma<<endl;
// cut on photopeak, +/- n*sigma
// TString cut = Form("data.LYSOChn>%f && data.LYSOChn<%f",
// peak_mean-sigma_cut*peak_sigma,
// peak_mean+sigma_cut*peak_sigma );
//cut off LYSO pedestal only
float lyso_pedestal_cut = 150;
// TString cut = Form("data.LYSOChn>%f ",lyso_pedestal_cut);//peak_mean-peak_sigma/2.0-100)
TString cut = Form("LYSOChn>%f ",lyso_pedestal_cut);//peak_mean-peak_sigma/2.0-100)
// TString cut2 = Form(" && data.LYSOChn<%f",peak_mean+peak_sigma/2.0-4.0*(peak_mean-lyso_pedestal_cut)/5.0);
TString cut2 = Form(" && LYSOChn<%f",peak_mean+peak_sigma/2.0-4.0*(peak_mean-lyso_pedestal_cut)/5.0);
// TString cut3 = Form(" && data.QuartzChn>%f",quartz_ped_cut);
TString cut3 = Form(" && QuartzChn>%f",quartz_ped_cut);
cout<<"cut : "<<cut<<endl;
cout<<"cut2: "<<cut2<<endl;
c1->cd(3);
// dTree->Draw("data.LYSOChn >> h4",cut + cut2);
dTree->Draw("LYSOChn>>h_lyso","");
TH1F* lyso_data = (TH1F*)h4->Clone("lyso_data");
lyso_data->Draw();
eff.lyso_hits = lyso_data->Integral(150,4096);
cout<<"Integral lyso data: "<<eff.lyso_hits<<endl;
c1->cd(2);
// dTree->Draw("data.QuartzChn>>h_quartz(4096,0,4096)","");
dTree->Draw("QuartzChn>>h_quartz(4096,0,4096)","");
c1->cd(4);
// dTree->Draw("data.QuartzChn >> h3",cut + cut2 + cut3);
dTree->Draw("QuartzChn >> h3",cut + cut2 + cut3);
TH1F* quartz_data = (TH1F*)h3->Clone("quartz_data");
quartz_data->Draw();
eff.quartz_hits = quartz_data->Integral(quartz_ped_cut,4096);
cout<<"Integral quartz data: "<<eff.quartz_hits<<endl;
eff.hit_efficiency = (float)eff.quartz_hits / eff.lyso_hits;
cout<<"Hit efficiency: "<< eff.hit_efficiency<<endl;
////////////////////////////////////////////
// PMT521, 10.44 ch/PE, 2.294pC/PE
// pedestal mean from run10965: 160.8+/-0.9, pedestal right edge: 165
// pedestal from beam-on runs: 165.2+/-1.0
// PE spectrum range: [0,377], where (4096-160.8)/10.44->377
// dTree->Draw("(data.QuartzChn-165.2)/10.44>>PE_h1(4000,0,377)",cut2);
TH1D *PE_h1 = new TH1D("PE_h1","",10000,0,(4096-quartz_pedestal)/OnePE);
TString PE_Cal = Form("(data.QuartzChn-%f)/%f>>PE_h1",quartz_pedestal,OnePE);
TString PE_Cal = Form("(QuartzChn-%f)/%f>>PE_h1",quartz_pedestal,OnePE);
dTree->Draw(PE_Cal, cut + cut3);
eff.total_PE = PE_h1->GetMean()*PE_h1->Integral();
eff.PE_eff = eff.total_PE/eff.lyso_hits;
// Light yield error: for y = x/a with OnePE = a+/-da = 10.44+/-sqrt(10.44), Delta y = (x/a^2)*da = x*a^(-3/2)
eff.PE_eff_err = eff.PE_eff/OnePE*sqrt(OnePE);
cout<<"====>> PE_h1->Integral(): "<<PE_h1->Integral()<<", PE_h1->GetMean():"<<PE_h1->GetMean()<<endl;
cout<<"====>> Total PE: "<<eff.total_PE<<", PE/hit:"<<eff.PE_eff<<" +/- "<< eff.PE_eff_err <<endl;
// delete h_lyso;
// delete s;
// delete fitF;
// delete h4;
// delete h3;
// delete PE_h1;
f_data->Close();
c1->SaveAs("Single1test.jpg");
}
what could the issue be?
Thank you for reading.