Efficiency function

Mansi · August 8, 2011, 8:00am

Hello

I have defined a cubic efficiency function in my code. And im fitting it with known efficiency. Then im evaluating efficiency using the same function at all values in a pre-defined histogram ( which is basically a energy spectrum).And then im dividing my original energy spectrum / histogram with those values. For that purpose i wrote the code given below:

[color=#0080BF]if(gamma_energy_KeV[i][j]>25.0)
{// cout<<"i "<<i;

   string fName;   
     sprintf(filename,"eff_parm_%i_%i.dat",i,j);
    fName = filename;
  
  infile.open(fName.c_str());           
   cout<<"File opened  "<<filename<<endl;

while(!infile.eof()||infile.good())
{infile>>a_1>>a_2>>a_3>>a_4;
fun->SetParameter(0,a_1);
fun->SetParameter(1,a_2);
fun->SetParameter(2,a_3);
fun->SetParameter(3,a_4);
// cout<<“a1:”<<a_1<<“a2:”<<a_2<<“a3:”<<a_3<<“a4:”<<a_4<<endl;
} infile.close();
value_eff[i][j]=fun->Eval(gamma_energy_KeV[i][j]); [color=#0040BF]// calculating the value of efficiency[/color]

[color=#0080FF]for (int d=0;d<4000;d++) // i have defined my histograms with a range of 0 to 4000
{ cout<<"i “<<i<<” j "<<j<<endl;
bin_count[i][j] =h_addback_Ge[i][j]->GetBinContent(d+1);
cout<<bin_count[i][j]<<“bin”<<endl;
corrected_bin_count[i][j]=bin_count[i][j]/value_eff[i][j];
h_addback_eff[i][j]->SetBinContent(d+1,corrected_bin_count[i][j]);
}[/color]

But the following code is showing a error. It does not calculate a correct value for the bin_count and hence also the corrected bin counts.

Kindly Help

Mansi

moneta · August 9, 2011, 3:42pm

Hi,
Can you please clarify better what the problem you are having is and eventually post the simplest possible running code with the required data file to run it

Best Regards

Lorenzo

Mansi · August 10, 2011, 5:32am

Hi

This was actually part of a code of the entire code. There are differents parts to the above code:

I have created a cubic function to fit the efficiency values of a detector with respect to energy.
The defined function for it is as follows:

[color=#8040FF]Double_t fun(Double_t *a,Double_t *par)
{
Double_t value1=0.0;
if (a[0]!=0)
{

value1=par[0]+(par[1]/a[0])+(par[2]/(a[0]*a[0]))+(par[3]/(a[0]*a[0]*a[0]));
return value1;

}}[/color]

Next in my code i have some predefined histograms which are basically energy spectrum . They are defined as:
[color=#8040FF]for(int i=0;i<5;i++)
{ for (int j=0;j<7;j++)
{
corrected_bin_count[i][j]=0.0;
bin_count[i][j]=0.0;
sprintf(name,“Ge_cluster%icrystal%i”,i,j);
h_GeEnergy[i][j] = new TH1F(name,name,4000.0,0.,4000.);
}}[/color]
Next I have to create a efficiency corrected energy spectrum . So for that purpose i wrote the following snippet:

First i get the value for efficiency at all energy points in the histograms by using the Eval () function.
[color=#8040FF]string fName;
sprintf(filename,“eff_parm_%i_%i.dat”,i,j);
fName = filename;
infile.open(fName.c_str());
cout<<"File opened "<<filename<<endl;
while(!infile.eof()||infile.good())
{infile>>a_1>>a_2>>a_3>>a_4;
fun->SetParameter(0,a_1);
fun->SetParameter(1,a_2);
fun->SetParameter(2,a_3);
fun->SetParameter(3,a_4);
// cout<<“a1:”<<a_1<<“a2:”<<a_2<<“a3:”<<a_3<<“a4:”<<a_4<<endl;
}
infile.close();

value_eff[i][j]=fun->Eval(gamma_energy_KeV[i][j]);
[/color]

Next i create a eff corrected spectrum , by dividing the counts by eff_value .
// there is another for loop over this which runs over i and j
for (int d=0;d<4000;d++)
{ bin_count[i][j] =h_GeEnergy[i][j]->GetBinContent(d+1);
cout<<bin_count[i][j]<<“bin”<<endl;
corrected_bin_count[i][j]=bin_count[i][j]/value_eff[i][j];
h_eff_corrected[i][j]->SetBinContent(d+1,corrected_bin_count[i][j]);
}

BUT doing such steps is nt giving desired results . It gives me 0 value for the bin_count[i][j].

Please Help…

moneta · August 10, 2011, 9:21am

Hi,
The code seems fine, but what is zero, is h_GeEnergy[i][j]->GetBinContent(d+1) ?
You should check how you have filled the h_GeEnergy histogram

Lorenzo

Mansi · August 10, 2011, 10:33am

Yes it gives me zero value for the GetBinContent . Yes i have properly filled the Ge Energy histogram…

moneta · August 10, 2011, 11:04am

If this is the case, you should check carefully that histograms, if the bin content is zero mens that that bin has not been filled with any entry

Lorenzo

Mansi · August 17, 2011, 9:09am

Hello

I have looked up into my code again… I seem to get stuck at a point. My question is now as follows:

I have originally defined a tree in which i have all my events , which im filling up in subsequent histograms in the program. Now should i first fill all the histograms and then calculate the efficiency value or can i do it simultaneouly.

moneta · August 19, 2011, 9:27pm

Hi,
I am not sure I have understood your problem. First you need to fill the histogram from the tree and then calculate your efficiency, but you can do this in a same application.

Lorenzo

Mansi · August 20, 2011, 8:23am

Hi…

I wanted to ask that , is it possible that to simultaneously fill a histogram event by event and also use that histogram to calculate the efficiency.

Mansi · August 21, 2011, 3:15am

Hello

I have attached a part of my entire code. Please look at the logic now , and tell me if its correct.

The programs breaks when it enters the loop to find the corrected bin count and fill a new histogram called h_addback_eff.

[color=#404040]The error which it shows is :::

*** Break *** segmentation violation
(no debugging symbols found)
Attaching to program: /proc/4757/exe, process 4757
(no debugging symbols found)…done.
done.
done.
done.
done.
[Thread debugging using libthread_db enabled]
[New Thread 0xb80556d0 (LWP 4757)]
done.
done.
done.
done.
done.

0x00110416 in __kernel_vsyscall ()
Missing separate debuginfos, use: debuginfo-install expat-2.0.1-5.i386 fontconfig-2.6.0-3.fc10.i386 freetype-2.3.7-1.fc10.i386 libX11-1.1.4-5.fc10.i386 libXau-1.0.4-1.fc10.i386 libXcursor-1.1.9-3.fc10.i386 libXdmcp-1.0.2-6.fc10.i386 libXext-1.0.4-1.fc9.i386 libXfixes-4.0.3-4.fc10.i386 libXft-2.1.13-1.fc10.i386 libXpm-3.5.7-4.fc9.i386 libXrender-0.9.4-3.fc9.i386 libjpeg-6b-43.fc10.i386 libpng-1.2.31-2.fc10.i386 libxcb-1.1.91-5.fc10.i386 zlib-1.2.3-18.fc9.i386
#1 0x003ecf43 in __waitpid_nocancel () from /lib/libc.so.6
#2 0x0038732b in do_system (
at …/sysdeps/posix/system.c:149
#3 0x0034259d in system (
at pt-system.c:29
#4 0x00876d8d in TUnixSystem::Exec ()
from /home/mansi/Download/root/lib/libCore.so
#5 0x0087c2ab in TUnixSystem::StackTrace ()
from /home/mansi/Download/root/lib/libCore.so
#6 0x0087d04d in TUnixSystem::DispatchSignals ()
from /home/mansi/Download/root/lib/libCore.so
#7 0x0087d14d in SigHandler () from /home/mansi/Download/root/lib/libCore.so
#8 0x008739f2 in sighandler () from /home/mansi/Download/root/lib/libCore.so
#9
#10 0x093adc55 in G__G__Hist_96_0_42 ()
from /home/mansi/Download/root/lib/libHist.so
#11 0x00e17d26 in Cint::G__ExceptionWrapper ()
from /home/mansi/Download/root/lib/libCint.so
#12 0x00e2b14f in G__exec_asm () from /home/mansi/Download/root/lib/libCint.so
#13 0x00f00946 in G__exec_loop () from /home/mansi/Download/root/lib/libCint.so
#14 0x00efe353 in G__exec_statement ()
from /home/mansi/Download/root/lib/libCint.so
#15 0x00efb15a in G__exec_statement ()
from /home/mansi/Download/root/lib/libCint.so
#16 0x00f002ef in G__exec_loop () from /home/mansi/Download/root/lib/libCint.so
#17 0x00efe353 in G__exec_statement ()
from /home/mansi/Download/root/lib/libCint.so
#18 0x00f002ef in G__exec_loop () from /home/mansi/Download/root/lib/libCint.so
#19 0x00efe353 in G__exec_statement ()
from /home/mansi/Download/root/lib/libCint.so
#20 0x00f002ef in G__exec_loop () from /home/mansi/Download/root/lib/libCint.so
#21 0x00efe353 in G__exec_statement ()
from /home/mansi/Download/root/lib/libCint.so
#22 0x00ea67ed in G__interpret_func ()
from /home/mansi/Download/root/lib/libCint.so
#23 0x00e946f7 in G__getfunction ()
from /home/mansi/Download/root/lib/libCint.so
#24 0x00e66f28 in G__getitem () from /home/mansi/Download/root/lib/libCint.so
#25 0x00e6ca73 in G__getexpr () from /home/mansi/Download/root/lib/libCint.so
#26 0x00efc82a in G__exec_statement ()
from /home/mansi/Download/root/lib/libCint.so
#27 0x00e523d2 in G__exec_tempfile_core ()
from /home/mansi/Download/root/lib/libCint.so
#28 0x00e526a9 in G__exec_tempfile_fp ()
from /home/mansi/Download/root/lib/libCint.so
#29 0x00f050b0 in G__process_cmd ()
from /home/mansi/Download/root/lib/libCint.so
#30 0x00866304 in TCint::ProcessLine ()
from /home/mansi/Download/root/lib/libCore.so
#31 0x00784e37 in TApplication::ProcessLine ()
from /home/mansi/Download/root/lib/libCore.so
#32 0x002e478e in TRint::HandleTermInput ()
from /home/mansi/Download/root/lib/libRint.so
#33 0x002e4245 in TTermInputHandler::Notify ()
from /home/mansi/Download/root/lib/libRint.so
#34 0x002e6a24 in TTermInputHandler::ReadNotify ()
from /home/mansi/Download/root/lib/libRint.so
#35 0x0087a439 in TUnixSystem::CheckDescriptors ()
from /home/mansi/Download/root/lib/libCore.so
#36 0x0087a799 in TUnixSystem::DispatchOneEvent ()
from /home/mansi/Download/root/lib/libCore.so
#37 0x007ebdd1 in TSystem::InnerLoop ()
from /home/mansi/Download/root/lib/libCore.so
#38 0x007eeb9b in TSystem::Run () from /home/mansi/Download/root/lib/libCore.so
#39 0x007828a7 in TApplication::Run ()
from /home/mansi/Download/root/lib/libCore.so
#40 0x002e6504 in TRint::Run () from /home/mansi/Download/root/lib/libRint.so
#41 0x08048db5 in main ()
The program is running. Quit anyway (and detach it)? (y or n) [answered Y; input not from terminal]
Detaching from program: /proc/4757/exe, process 4757
Root > Function histograms() busy flag cleared[/color]

[color=#004080]gROOT->Reset();
#include
#include
#include <TFile.h>
#include <TTree.h>
#include <TCutG.h>
#include
#include <TROOT.h>
#include
#include “snprintf.h”

using namespace std;

void histograms()
{

Float_t SiEnergy[6];
Float_t SiTime[6];
Float_t GeEnergy[5][7];
Float_t GeTime[5][7] ;
Float_t MagFieldUp;
Float_t MagFieldDown;
Float_t gamma_energy_KeV[5][7];
Float_t gamma_DGF_time[5][7];
Float_t silicon_energy[6];
Float_t silicon_time[6];
Float_t gamma_doppler[5][7];
Float_t ge_energy[5][7];
Float_t mag_up;
Float_t mag_down;
Float_t Ge_mult[5][7];
Float_t GeEnergy_1[5][7];//these have energy filled in for mult =2
Float_t GeEnergy_2[5][7];
Int_t gamma_mult;
Int_t crystal_hit[5][7];
Int_t x,z,x2,z2;
Float_t phi_ge,phi_si,test;
Float_t corrected_counts[5][7],corrected_bin_count[5][7],bin_count[5][7];
int bin_num;
char name[100];
TH1F *h_SiEnergy[6];
TH1F *h_GeEnergy[5][7];
TH1F *h_addback_Ge[5][7];
TH1F *h_addback_eff[5][7];
TH1F *h_addback_cluster_Ge[5];
TH1F *h_gamma_gated[5][7];
TH1F *h_gamma1;
TH1F *h_SiTime[6];
TH1F *h_gamma2;
TH2F *h_clock_time;
TH2F *h_Si_Ge_addback;
TH2F *h_Ge_Ge;
TH1F *h_coin_gamma;
TH1F *h_add_all;
ifstream infile;
Double_t value_eff[5][7];

TF1 *fun = new TF1(“fit”,fit,210.0,1600,4);

for(int i=0;i<5;i++)
{ for (int j=0;j<7;j++)
{gamma_doppler[i][j]=0.0;
corrected_bin_count[i][j]=0.0;
bin_count[i][j]=0.0;
sprintf(name,“Ge_cluster%icrystal%i”,i,j);
h_GeEnergy[i][j] = new TH1F(name,name,4000.0,0.,4000.);
}}

for(int i=0;i<5;i++)
{ for (int j=0;j<7;j++)
{sprintf(name,“Ge_cluster%icrystal%i_gated”,i,j);
h_gamma_gated[i][j] = new TH1F(name,name,4000.0,0.,4000.);
}}

for(int i=0;i<5;i++)
{
for (int j=0;j<7;j++)
{ sprintf(name,“Ge_addback_cluster%icrystal%i”,i,j);
//sprintf(name,“h_Si_en_Ge_time[%i]”,i);
h_addback_Ge[i][j] = new TH1F(name,name,4000.0,0.,4000.);
}}

TFile *f2=new TFile(“FinalOutput_de_tdc.root”);
TTree miT =(TTree)f2->Get(“miT”);
// TTree T1 = (TTree)f1->Get(“T1”);
miT->SetBranchAddress(“gamma_energy_KeV”,&gamma_energy_KeV);
miT->SetBranchAddress(“gamma_DGF_time”,&gamma_DGF_time);
miT->SetBranchAddress(“silicon_energy”,&silicon_energy);
miT->SetBranchAddress(“silicon_time”,&silicon_time);

for(Int_t l=0;l<5;l++)
{ for(Int_t n=0;n<7;n++)
{ crystal_hit[l][n] = 0;
ge_energy[i][j]=0.0;
}}
Int_t n_entries=(Int_t)miT->GetEntries();
cout<<“n_entries =”<<n_entries<<endl;
cout<<“Entering loop”<<endl;
float phi;

//STARTING OF THE LOOP ON ENTRY IN TREE
//*************************************
for (int k=0;k<n_entries;k++)
{//cout<<"**"<<endl;
miT->GetEntry(k);
// T1->GetEntry(k);
gamma_mult =0;

char filename[200];

for(int i=0;i<5;i++)
{ for( int j=0;j<7;j++){

  h_addback_Ge[i][j]->Fill(gamma_energy_KeV[i][j]);
       h_add_all->Fill(gamma_energy_KeV[i][j]);
    
   //EFFICIENCY CALCULATION BEGINS HERE

if(gamma_energy_KeV[i][j]>25.0)
{ string fName;
sprintf(filename,“eff_parm_%i_%i.dat”,i,j);
fName = filename;
infile.open(fName.c_str());
cout<<“File opened “<<filename<<endl;
while(!infile.eof()||infile.good())
{infile>>a_1>>a_2>>a_3>>a_4;
fun->SetParameter(0,a_1);
fun->SetParameter(1,a_2);
fun->SetParameter(2,a_3);
fun->SetParameter(3,a_4);
// cout<<“a1:”<<a_1<<“a2:”<<a_2<<“a3:”<<a_3<<“a4:”<<a_4<<endl;
}
infile.close();
value_eff[i][j]=fun->Eval(gamma_energy_KeV[i][j]);
cout<<value_eff[i][j]<<”<<<”<<i<<j<<endl;
cout<<gamma_energy_KeV[i][j]<<endl;
corrected_counts[i][j]=gamma_energy_KeV[i][j]/value_eff[i][j];

for (int d=0;d<4000;d++)
{ bin_count[i][j] =h_addback_Ge[i][j]->GetBinContent(d+1);
if (bin_count[i][j])>0.0{
cout<<bin_count[i][j]<<" bin"<<endl;
corrected_bin_count[i][j]=bin_count[i][j]/value_eff[i][j];
h_addback_eff[i][j]->SetBinContent(d+1,corrected_bin_count[i][j]);}
}} //for loop and if loop

}}

c3->cd();
h_addback_eff[0][0]->Draw();
c4->cd();
h_addback_Ge[0][0]->Draw();
c5->cd();
h_addback_eff[0][1]->Draw();
c6->cd();
h_addback_Ge[0][1]->Draw();
c7->cd();
h_addback_eff[0][2]->Draw();
c8->cd();
h_addback_Ge[0][2]->Draw();
c9->cd();
h_addback_eff[0][3]->Draw();

}

Double_t fit(Double_t *a,Double_t *par)
{
Double_t value1=0.0;
Double_t a_ene=0.0,b_ene=0.0;
// cout<<“ee”<<endl;
if (a[0]!=0)
{
value1=par[0]+(par[1]/a[0])+(par[2]/(a[0]*a[0]))+(par[3]/(a[0]*a[0]*a[0]));
return value1;
}}

[/color]

Mansi · August 24, 2011, 7:14am

Please look into the error and the code which i have attached…

moneta · August 24, 2011, 10:16am

Please upload the code as attachment together with the required files needed to run your code, so we can have a look into it

Lorenzo

Mansi · August 24, 2011, 12:18pm

In my code i have first build a Tree and stored all the events in it , after doing some initial calibration. Then in a different code ( where iam calculating efficiency ) i use that Tree to get the events. Its not possible to attach tht tree ( its pretty huge file .) .

Suggest some other way !!!

moneta · August 24, 2011, 3:04pm

By looking at the code, should not get the efficiency corrected spectrum by filling the corrected histogram using the efficiency as a weight ?
So, basically doing

histo->Fill(gamma_energy, efficiency )

From your code it seems you are mixing the x quantities ( the energy) with the bin contents (the counts /bin)

Lorenzo

Mansi · August 24, 2011, 10:07pm

Yes , i guess you are right … Can you suggest a optimum way to do it …