ROOT macro for calculation of neutron rates in the BRIKEN detector

Here a ROOT macro to calculate neutron rates from onlines root files generated DAQ configuration 
160718Conf_BrikenFull_cal5.xlsx (and the later ones).

Setup parameters are indicated at the begining of the macro. An output data file is generated with the rates by counter, ring and ADC. If the input root file is "file.root", the 
output data file is "file.root_data".

Further questions by email to: atarisal@gmail.com


 

//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();

}