AnalyseCombinedGammaFiles.C

Go to the documentation of this file.

// created on 23.09.13 by Marcell Steinen steinen@kph.uni-mainz.de
// root macro to analyse the combined energy spectrum
//full chain description:
// root macro to create an energy spectrum of parts of simulations with a high number of events (memory leak somewhere :/); to be used with "script_Gamma_Ana_List_split.sh" on himster
// further analysis is done  via bash script "script_CombineAndAnalyzeGamma.sh" which combines the split up root output via "hadd" and calls root macro "AnalyseCombinedGammaFiles.C" to analyse the spectrum

#include<TH1D.h>
#include<TH2D.h>
#include<stdio.h>
#include<TString.h>
#include"TFile.h"
#include<iostream>
#include<fstream>
#include"TROOT.h"

using namespace std;

Double_t Compton(Double_t Energy, Double_t Th)
{
return Energy-Energy/(1+Energy/0.000511*(1-TMath::Cos(Th)));
}
Double_t invCompton(Double_t Energy, Double_t Eprime)
{
return TMath::ACos((Energy*Energy-Energy*Eprime-Eprime*0.000511)/(Energy*(Energy-Eprime)));
}
Double_t PeakFunc (Double_t *x, Double_t *par)
{
        return par[0]+x[0] * par[1]+x[0] * x[0] * par[2]+ par[3] * 1/(TMath::Sqrt(2 * TMath::Pi()) * par[5]) * TMath::Exp(- pow((x[0] -par[4]),2)/2/par[5]/par[5]); 
}
Double_t PoissonFunc(Double_t *x, Double_t *par)
{
        return par[1]*TMath::Poisson(x[0],par[0]);
}

int AnalyseCombinedGammaFiles(TString Filepath, TString FileName,Int_t NoOfJobs = 10)//TString FileListname_ext,TString OutfileName_ext)
{
        
        TString Anapath = Filepath;
        Anapath.Remove(Anapath.Last('/'),Anapath.Length()-Anapath.Last('/'));
        cout << Anapath << endl;
        TString InputName =Filepath.Remove(0,Filepath.Last('/')+1);
        
        if (InputName.Contains("TripleBall40Offset10STT"))
                TString Geometry = "35";
        if (InputName.Contains("TripleBall40Offset20STT"))
                TString Geometry = "36";        
        if (InputName.Contains("SecTar"))
                Int_t Sektar = 1;
        else
                Int_t Sektar = 0;
        
        TFile *File = new TFile(FileName.Data());
        File->ls();
        TH1D *hGammaSpectrum;
        TH1D *hNoOfHits;
        TH1D *hParBuffer;
        File->GetObject("gamTde;1",hGammaSpectrum);
        File->GetObject("Number of Hits;1",hNoOfHits);
        File->GetObject("ParameterBuffer;1",hParBuffer);
        Double_t Energy = hParBuffer->GetBinContent(1)/NoOfJobs;
        Int_t nEvents = hParBuffer->GetBinContent(2)/NoOfJobs;
        Int_t ChannelResolution = hParBuffer->GetBinContent(3)/NoOfJobs;
        Double_t Resolution = hParBuffer->GetBinContent(4)/NoOfJobs;
        cout << Energy << "  " << nEvents << "   " << ChannelResolution << "   " << Resolution << endl;
        
        ofstream txtfile;
        txtfile.open("OutputCombinedAna.txt");
        txtfile << "File analysed:\t" << InputName << endl;
        txtfile<< "Number of Simulated Events:\t"<<nEvents<<endl;
        
        
        
        TCanvas* can3 = new TCanvas("can3","germanium detector",0,0,1000,1000);
        hGammaSpectrum->Draw();
  
  
  //Analysis of spectrum
  
        Double_t lowAngle = invCompton(0.001332,0.001040);                      //get the angle from the values of co60 for comparison
        Double_t highAngle = invCompton(0.001332,0.001096);                     //get the angle from the values of co60 for comparison
        Int_t npeaks = 1;
        Int_t PeakToLook = 1;
        
        TSpectrum *s = new TSpectrum(npeaks);
        s->Search(hGammaSpectrum,npeaks,"new",0.01);

        Float_t *xpeaks = s->GetPositionX();
        Float_t *ypeaks = s->GetPositionY();
        
                //Print Peaks
        for (Int_t i = 0; i < npeaks; ++i)
        { 
                cout <<"("<< xpeaks[i]<<"," << ypeaks[i] << ") ";
        }
        cout << endl;
        
        Double_t PeakX = xpeaks[PeakToLook-1]*ChannelResolution;
        Double_t PeakY = ypeaks[PeakToLook-1];
        Double_t DPeakY = sqrt(PeakY); 
        Double_t SumPeak = 0;
        
        for(Int_t i = PeakX-Resolution*10*10/2.3548; i <= PeakX+Resolution*10*10/2.3548; i++)           //Peak = 10 * Resolution of 1 crystal --> takes multiple hits into account (faster than fitting and than doing it again, error is small)
                SumPeak += hGammaSpectrum->GetBinContent(i);
        Double_t DSumPeak = sqrt(SumPeak);
        cout <<"MaxX: "<< PeakX << "; MaxY: " << PeakY << " +- " << DPeakY<< endl;
        cout << "SumPeak: " << SumPeak << " +- " << DSumPeak << endl;
        txtfile <<"MaxX: "<< PeakX << "; MaxY: " << PeakY << " +- " << DPeakY<< endl;
        txtfile << "SumPeak: " << SumPeak << " +- " << DSumPeak << endl;
        
        Double_t ComptonEdge = Int_t(PeakX*(1-1/(1+2*PeakX/5110)))+1;
        cout <<"ComptonEdge @ Bin: " << ComptonEdge << endl;
        Double_t SumCompton=0, DSumCompton = 0; 
        
        
  for (Int_t i =Compton(Energy,lowAngle)*ChannelResolution; i <= Compton(Energy,highAngle)*ChannelResolution; i++)
  {
                SumCompton += hGammaSpectrum->GetBinContent(i);
        }
        DSumCompton = sqrt(SumCompton);
        
        cout << "lowEnChan: " << Compton(Energy,lowAngle)*ChannelResolution << endl << "highEnChan: "<< Compton(Energy,highAngle)*ChannelResolution << endl;
        txtfile << "lowEnChan: " << Compton(Energy,lowAngle)*ChannelResolution << endl << "highEnChan: "<< Compton(Energy,highAngle)*ChannelResolution << endl;
        cout <<"N. of Compton Events: "<< SumCompton<< " +- " << DSumCompton << endl;
        txtfile <<"N. of Compton Events: "<< SumCompton<< " +- " << DSumCompton << endl;
        Double_t PeakToCompton = 0;
        Double_t DPeakToCompton = 0;
        if (SumCompton != 0)
        {
                PeakToCompton = PeakY/SumCompton*(Compton(Energy,highAngle)-Compton(Energy,lowAngle))*ChannelResolution;
                DPeakToCompton = sqrt(pow(DPeakY/SumCompton,2)+ pow(PeakY*DSumCompton/SumCompton/SumCompton,2))*(Compton(Energy,highAngle)-Compton(Energy,lowAngle))*ChannelResolution;
                cout << "Peak/Compton (res): "<< PeakToCompton << " +- "<< DPeakToCompton << endl; 
                txtfile << "Peak/Compton (res):\t"<< PeakToCompton << " +- "<< DPeakToCompton << endl;
        }
        else
        {
                
                cout << "SumCompton = 0 -> no PeakToCompton possible!"<<endl;
                txtfile << "SumCompton = 0 -> no PeakToCompton possible!"<<endl;
        }
        Double_t AllEntries = hGammaSpectrum->GetEntries();
        cout << "Entries: " << AllEntries << endl;
        txtfile << "Entries: " << AllEntries << endl;
        
        
        //fitting the histograms
        
        TF1 *GausBG = new TF1("GausBG","gausn"/*(0)+pol1(3)" */,(PeakX-100)/ChannelResolution,(PeakX+100)/ChannelResolution);
                GausBG->SetParName(0,"Ampl");
                GausBG->SetParName(1,"x0");
                GausBG->SetParName(2,"sigma");
                GausBG->SetParameter(0,PeakY/10000);
                GausBG->SetParameter(1,PeakX/Double_t(ChannelResolution));
                GausBG->SetParameter(2,100./Double_t(ChannelResolution));
                GausBG->SetParLimits(0,0.0000001, 1000000);
                GausBG->SetParLimits(1,(PeakX-10)/Double_t(ChannelResolution),(PeakX+10)/Double_t(ChannelResolution));
                GausBG->SetParLimits(2,1./Double_t(ChannelResolution),100./Double_t(ChannelResolution));
                /*GausBG->SetParameter(0,1);
                GausBG->SetParameter(1,1);
                GausBG->SetParameter(2,0);
                GausBG->FixParameter(2,0);
                GausBG->FixParameter(1,0);
                GausBG->FixParameter(0,0);
                */
                GausBG->SetNpx(100000);
                GausBG->SetLineColor(kRed);     
        hGammaSpectrum->Fit(GausBG,"R");
        
        txtfile << "Fitparameter Spectrum (model: gausn):" << endl;
        for (Int_t i = 0; i < 3; i++)
        {
                txtfile << GausBG->GetParName(i) << ":\t" << GausBG->GetParameter(i) << endl;
        }
        Double_t FWHM = GausBG->GetParameter(2)*2.3548200;
        txtfile << "FWHM[keV]:\t" <<  FWHM*1000000 << endl;
        cout << "FWHM[keV]:\t" <<  FWHM*1000000 << endl;
        
        /*TF1 *Poisson = new TF1("Poisson","gausn",0,15);//,2);                 //maybe gausn? maybe no fit? //gausn --> no ,2);
                //Poisson->SetParameter(0,4);
                //Poisson->SetParameter(1,SumPeak);
                //Poisson->SetParLimits(0,0,12);
                Poisson->SetParName(0,"Ampl");//Poisson->SetParName(0,"lambda");
                Poisson->SetParName(1,"x0");//Poisson->SetParName(1,"Ampl");
                Poisson->SetParName(2,"sigma");
                //Poisson->SetParameter(1,1000);
                Poisson->SetNpx(100000);
                Poisson->SetLineColor(kRed);    
        //hNoHits->Fit(Poisson);*/
        
        //txtfile << "Fitparameter Number of Hits (model: Ampl.*Poisson):" << endl;
        //for (Int_t i = 0; i < 3; i++)                 // if gausn fit --> < 2 -> < 3 , if no fit, comment it out!!!
        //{
        //      txtfile << Poisson->GetParName(i) << ":\t" << Poisson->GetParameter(i) << endl;
        //}
        
        //some make-up for the histograms
        Double_t FullEnergyPeakEff = double(int(GausBG->GetParameter(0)*ChannelResolution/(nEvents*2)*100*1000))/1000;
        Double_t FullEnergyPeakEffError = double(int(GausBG->GetParError(0)*ChannelResolution/(nEvents*2)*100*1000))/1000;
        
        cout << "Full-Energy-Peak-Eff. [%]: " << FullEnergyPeakEff << endl;; // *2 because only 2Pi of solid angle is simulated
        txtfile << "Full-Energy-Peak-Eff. [%]:\t" << FullEnergyPeakEff << endl; // *2 because only 2Pi of solid angle is simulated 
        // writing to files and closing
        cout << "Error ofFull-Energy-Peak-Eff. [%]: " <<  FullEnergyPeakEffError<< endl; // *2 because only 2Pi of solid angle is simulated
        txtfile << "Error of Full-Energy-Peak-Eff. [%]:\t" << FullEnergyPeakEffError << endl; // *2 because only 2Pi of solid angle is simulated 
        // writing to files and closing 
        TFile *Output = new TFile("OutputCombinedAna.root","RECREATE");
        hGammaSpectrum->Write();
        hNoOfHits->Write();
        
        txtfile.close();
        Output->Close();
        File->Close();
        
        TString CollectFileName = Anapath +"/CollectFEPF.txt";
        ofstream CollectFile(CollectFileName, ios::out|ios::app);                       //collect name and Full Energy Eff of all Sims here
                CollectFile <<  Geometry << "\t" << Sektar <<  "\t"<< Energy << "\t"<< nEvents << "\t" << FullEnergyPeakEff << "\t" << FullEnergyPeakEffError << "\t" <<FWHM*1000000 << "\t"<< PeakToCompton<< "\t"<<InputName << endl;
        
        CollectFile.close();
        
        cout << "macro finished successfully" << endl;
        
  return 0;
}