//ROOT macro for calculation of neutron rates in the BRIKEN detector
// A. Tarifeno-Saldivia && A. Tolosa, Date: 12Oct16, update 12Oct16
TH1F* h_sum = NULL;
void BRIKEN_rates_v3(){
gStyle->SetOptStat("nei");
//**************** SETUP ****************
int Ndet_i=1;
int Ndet_f=140;
Double_t deltaPulser=30.0;
Double_t PulserLowLimit=1350.0;
Double_t Emin = 160.0; // in keV
Double_t Emax = 900.0; // in keV
Double_t PeakAmplitude_threshold=0.005; // 0.005 should work.
// Higher values than 0.08 may affect finding pulser position
//*****************************************
TString input;//="0cm_2.root";
cout<<"Existing root files in current directory: \n";
gSystem->Exec("ls *.root");
cout<<"Input root file: \n";
cin >> input;
TFile * file = TFile::Open(input);
ofstream outFile(input + "_data");
TH1F * inter = (TH1F*) file->Get("He001");
h_sum = (TH1F*) inter->Clone();
h_sum -> Reset();
h_sum->SetName("Sum histogram");
h_sum->SetTitle("BRIKEN full histogram: " + input);
TString histName;
TCanvas *C1= new TCanvas("C1", "C1", 1);
C1->cd();
double Counter_nRates[141];
for(int i=0 ; i < 141 ; i++){Counter_nRates[i]=0.0;}
Double_t TotalnRate=0.0;
outFile<<"# File:\t"<<input<<endl;
outFile<<"# Histograms size:\t"<<h_sum->GetSize()<<endl;
outFile<<"# Emin:\t"<<Emin<<" keV\tEmax:\t"<<Emax<<" keV\n";
outFile<<"# deltaPulser:\t"<<deltaPulser<<" keV\tPulserLowLimit:\t"<<PulserLowLimit<<" keV\n";
outFile<<"#####################################################################################################\n";
outFile<<"# Detector\tnRate\tNeutronCounts\tPulserCounts\tpulserPos\tPeakBINo\tPeakBINf\n";
outFile<<"#\t\t(counts\t(counts)\t (counts)\t(keV)\t \t(bin)\t \t(bin)\t\n";
outFile<<"#\t\t/pulser)\n";
outFile<<"#####################################################################################################\n";
for(int i=Ndet_i; i<=Ndet_f ; i++){
if(i>0 && i<10) histName= Form("He00%d", i);
if(i>=10 && i<100) histName= Form("He0%d", i);
if(i>=100 && i<=140) histName= Form("He%d", i);
TH1F* hi =(TH1F*)file -> Get(histName);
// Find peak for pulser
TSpectrum *s = new TSpectrum(100);
Int_t nfound = s->Search(hi,10,"",PeakAmplitude_threshold);
TList *functions = hi->GetListOfFunctions();
TPolyMarker *pm = (TPolyMarker*)functions->FindObject("TPolyMarker");
Double_t* pmarrayX=pm->GetX() ; // Find polymarkers position
Double_t pulserPos, PulserCounts, NeutronCounts, nRate;
bool found=false;
int peakCounter=0;
int BinEmin= hi->GetXaxis()->FindBin(Emin);
int BinEmax= hi->GetXaxis()->FindBin(Emax);
for (int j = 0; j < nfound; j++) {
pulserPos=pmarrayX[j];
if(pulserPos>=PulserLowLimit){
peakCounter++;
found=true;
TF1 *f1 = new TF1("f1","gaus",-.5,.5);
f1->SetParameters(1000.0, pulserPos, deltaPulser);
hi -> Fit("f1", "Q", "",pulserPos-deltaPulser,pulserPos+deltaPulser);
Double_t mean = f1->GetParameter(1);
Double_t sd = f1->GetParameter(2);
Double_t Peak_o=mean-3.5*sd;
Double_t Peak_f=mean+3.5*sd;
int PeakBINo= hi->GetXaxis()->FindBin(Peak_o);
int PeakBINf= hi->GetXaxis()->FindBin(Peak_f);
PulserCounts = hi->Integral(PeakBINo,PeakBINf);
NeutronCounts = hi->Integral(BinEmin,BinEmax);
nRate=NeutronCounts/PulserCounts;
Counter_nRates[i]=nRate;
//TotalnRate+=nRate;
outFile<<histName<<"\t"<<nRate<<"\t"<<NeutronCounts<<"\t"<<PulserCounts<<"\t";
outFile<<pulserPos<<"\t"<<Peak_o<<"\t"<<Peak_f<<"\t"<<endl;
// cout<<histName<<"\t"<<nRate<<"\t"<<NeutronCounts<<"\t"<<PulserCounts<<"\t";
// cout<<pulserPos<<"\t"<<PeakBINo<<"\t"<<PeakBINf<<"\t"<<hi->GetSize()<<endl;
}
if(j==nfound-1 && found==false){
outFile<<histName<<"\t -------> ERROR FINDING THE PULSER, CHECK PeakAmplitude_threshold !!!!!\n";
cout<<histName<<"\t -------> ERROR FINDING THE PULSER, CHECK PeakAmplitude_threshold !!!!!\n";
}
if(j==nfound-1 && found==true && peakCounter>1){
outFile<<histName<<"\t -------> MORE THAN ONE PULSER POSITION, CHECK PeakAmplitude_threshold !!!!!\n";
cout<<histName<<"\t -------> MORE THAN ONE PULSER POSITION, CHECK PeakAmplitude_threshold !!!!!\n";
}
}
h_sum -> Add(hi);
}
outFile<<"\n#######################################################################\n";
outFile<<"# NEUTRON RATES (counts/pulser)\n";
outFile<<"#######################################################################\n";
//**************** Rates Calculation:
// ADC arrays
int V1A[6][16]={{1,2,3,4,5,6,7,8,9,10,11,12,21,22,23,24},
{13,14,15,16,17,18,19,20,31,32,33,34,35,36,37,38},
{25,26,27,28,29,30,43,44,45,46,47,48,49,50,51,52},
{39,40,41,42,53,54,55,56,57,58,59,60,61,62,63,64},
{65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80},
{81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,0}};
int V2A[6][8]= {{96,97,98,99,100,101,102,103},
{104,105,106,107,108,109,110,0},
{111,112,113,114,115,116,117,118},
{119,120,121,122,123,124,125,0},
{126,127,128,129,130,131,132,133},
{134,135,136,137,138,139,140,0}};
int Ri[7][28]={{1,2,3,7,8,9,21,13,14,15,19,20,26,27,28,29,30,46,47,48,49,50,0,0,0,0,0,0},
{4,5,10,11,22,23,16,17,33,34,35,36,53,54,55,56,0,0,0,0,0,0,0,0,0,0,0,0},
{6,12,24,18,31,32,37,38,25,43,44,45,51,52,39,40,41,42,57,58,59,60,61,62,63,64,0,0},
{66,67,74,75,81,82,83,84,85,111,112,113,114,115,0,0,0,0,0,0,0,0,0,0,0,0,0,0},
{65,68,73,76,86,87,88,89,90,91,116,117,118,119,120,121,0,0,0,0,0,0,0,0,0,0,0,0},
{92,93,94,95,96,97,98,99,100,122,123,124,125,126,127,128,129,130,0,0,0,0,0,0,0,0,0,0},
{69,70,71,72,77,78,79,80,101,102,103,104,105,106,107,108,109,110,131,132,133,134,135,136,137,138,139,140}};
//total
for(int i=0 ; i < 141 ; i++){TotalnRate+=Counter_nRates[i];}
cout<<"#########################################\n";
cout<<"# Total neutron rate (counts/pulser):\t"<<TotalnRate<<endl;
cout<<"#########################################\n";
cout<<" RATE BY RING:"<<endl;
cout<<"---------------"<<endl;
outFile<<"Total neutron rate:\t"<<TotalnRate<<endl;
outFile<<"#########################################\n";
outFile<<" BY RING:"<<endl;
outFile<<"#########################################\n";
for(int i=0 ; i < 7 ; i++){
double ringRATE=0.0;
for(int j=0 ; j < 28 ; j++){
ringRATE+=Counter_nRates[Ri[i][j]];
}
cout<<"R"<<i+1<<"\t"<<ringRATE<<endl;
outFile<<"R"<<i+1<<"\t"<<ringRATE<<endl;
}
outFile<<"#########################################\n";
outFile<<" BY ADC:"<<endl;
outFile<<"#########################################\n";
cout<<"#########################################\n";
cout<<" RATE BY ADC:"<<endl;
cout<<"---------------"<<endl;
//By ADC
//VME1
for(int i=0 ; i < 6 ; i++){
double vmeRATE=0.0;
for(int j=0 ; j < 16 ; j++){
vmeRATE+=Counter_nRates[V1A[i][j]];
}
cout<<"V1A"<<i+1<<"\t"<<vmeRATE<<endl;
outFile<<"V1A"<<i+1<<"\t"<<vmeRATE<<endl;
}
//VME2
for(int i=0 ; i < 6 ; i++){
double vmeRATE=0.0;
for(int j=0 ; j < 8 ; j++){
vmeRATE+=Counter_nRates[V2A[i][j]];
}
cout<<"V2A"<<i+1<<"\t"<<vmeRATE<<endl;
outFile<<"V1A"<<i+1<<"\t"<<vmeRATE<<endl;
}
h_sum -> Draw();
gPad->SetLogy();
outFile.close();
//gApplication->Terminate();
}
|