FairRoot/PandaRoot
Functions
AnalyseCombinedGammaFiles.C File Reference
#include <TH1D.h>
#include <TH2D.h>
#include <stdio.h>
#include <TString.h>
#include "TFile.h"
#include <iostream>
#include <fstream>
#include "TROOT.h"

Go to the source code of this file.

Functions

Double_t Compton (Double_t Energy, Double_t Th)
 
Double_t invCompton (Double_t Energy, Double_t Eprime)
 
Double_t PeakFunc (Double_t *x, Double_t *par)
 
Double_t PoissonFunc (Double_t *x, Double_t *par)
 
int AnalyseCombinedGammaFiles (TString Filepath, TString FileName, Int_t NoOfJobs=10)
 

Function Documentation

int AnalyseCombinedGammaFiles ( TString  Filepath,
TString  FileName,
Int_t  NoOfJobs = 10 
)

Definition at line 35 of file AnalyseCombinedGammaFiles.C.

References can3, Compton(), Double_t, i, invCompton(), nEvents, out, s, sqrt(), and TString.

36 {
37 
38  TString Anapath = Filepath;
39  Anapath.Remove(Anapath.Last('/'),Anapath.Length()-Anapath.Last('/'));
40  cout << Anapath << endl;
41  TString InputName =Filepath.Remove(0,Filepath.Last('/')+1);
42 
43  if (InputName.Contains("TripleBall40Offset10STT"))
44  TString Geometry = "35";
45  if (InputName.Contains("TripleBall40Offset20STT"))
46  TString Geometry = "36";
47  if (InputName.Contains("SecTar"))
48  Int_t Sektar = 1;
49  else
50  Int_t Sektar = 0;
51 
52  TFile *File = new TFile(FileName.Data());
53  File->ls();
54  TH1D *hGammaSpectrum;
55  TH1D *hNoOfHits;
56  TH1D *hParBuffer;
57  File->GetObject("gamTde;1",hGammaSpectrum);
58  File->GetObject("Number of Hits;1",hNoOfHits);
59  File->GetObject("ParameterBuffer;1",hParBuffer);
60  Double_t Energy = hParBuffer->GetBinContent(1)/NoOfJobs;
61  Int_t nEvents = hParBuffer->GetBinContent(2)/NoOfJobs;
62  Int_t ChannelResolution = hParBuffer->GetBinContent(3)/NoOfJobs;
63  Double_t Resolution = hParBuffer->GetBinContent(4)/NoOfJobs;
64  cout << Energy << " " << nEvents << " " << ChannelResolution << " " << Resolution << endl;
65 
66  ofstream txtfile;
67  txtfile.open("OutputCombinedAna.txt");
68  txtfile << "File analysed:\t" << InputName << endl;
69  txtfile<< "Number of Simulated Events:\t"<<nEvents<<endl;
70 
71 
72 
73  TCanvas* can3 = new TCanvas("can3","germanium detector",0,0,1000,1000);
74  hGammaSpectrum->Draw();
75 
76 
77  //Analysis of spectrum
78 
79  Double_t lowAngle = invCompton(0.001332,0.001040); //get the angle from the values of co60 for comparison
80  Double_t highAngle = invCompton(0.001332,0.001096); //get the angle from the values of co60 for comparison
81  Int_t npeaks = 1;
82  Int_t PeakToLook = 1;
83 
84  TSpectrum *s = new TSpectrum(npeaks);
85  s->Search(hGammaSpectrum,npeaks,"new",0.01);
86 
87  Float_t *xpeaks = s->GetPositionX();
88  Float_t *ypeaks = s->GetPositionY();
89 
90  //Print Peaks
91  for (Int_t i = 0; i < npeaks; ++i)
92  {
93  cout <<"("<< xpeaks[i]<<"," << ypeaks[i] << ") ";
94  }
95  cout << endl;
96 
97  Double_t PeakX = xpeaks[PeakToLook-1]*ChannelResolution;
98  Double_t PeakY = ypeaks[PeakToLook-1];
99  Double_t DPeakY = sqrt(PeakY);
100  Double_t SumPeak = 0;
101 
102  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)
103  SumPeak += hGammaSpectrum->GetBinContent(i);
104  Double_t DSumPeak = sqrt(SumPeak);
105  cout <<"MaxX: "<< PeakX << "; MaxY: " << PeakY << " +- " << DPeakY<< endl;
106  cout << "SumPeak: " << SumPeak << " +- " << DSumPeak << endl;
107  txtfile <<"MaxX: "<< PeakX << "; MaxY: " << PeakY << " +- " << DPeakY<< endl;
108  txtfile << "SumPeak: " << SumPeak << " +- " << DSumPeak << endl;
109 
110  Double_t ComptonEdge = Int_t(PeakX*(1-1/(1+2*PeakX/5110)))+1;
111  cout <<"ComptonEdge @ Bin: " << ComptonEdge << endl;
112  Double_t SumCompton=0, DSumCompton = 0;
113 
114 
115  for (Int_t i =Compton(Energy,lowAngle)*ChannelResolution; i <= Compton(Energy,highAngle)*ChannelResolution; i++)
116  {
117  SumCompton += hGammaSpectrum->GetBinContent(i);
118  }
119  DSumCompton = sqrt(SumCompton);
120 
121  cout << "lowEnChan: " << Compton(Energy,lowAngle)*ChannelResolution << endl << "highEnChan: "<< Compton(Energy,highAngle)*ChannelResolution << endl;
122  txtfile << "lowEnChan: " << Compton(Energy,lowAngle)*ChannelResolution << endl << "highEnChan: "<< Compton(Energy,highAngle)*ChannelResolution << endl;
123  cout <<"N. of Compton Events: "<< SumCompton<< " +- " << DSumCompton << endl;
124  txtfile <<"N. of Compton Events: "<< SumCompton<< " +- " << DSumCompton << endl;
125  Double_t PeakToCompton = 0;
126  Double_t DPeakToCompton = 0;
127  if (SumCompton != 0)
128  {
129  PeakToCompton = PeakY/SumCompton*(Compton(Energy,highAngle)-Compton(Energy,lowAngle))*ChannelResolution;
130  DPeakToCompton = sqrt(pow(DPeakY/SumCompton,2)+ pow(PeakY*DSumCompton/SumCompton/SumCompton,2))*(Compton(Energy,highAngle)-Compton(Energy,lowAngle))*ChannelResolution;
131  cout << "Peak/Compton (res): "<< PeakToCompton << " +- "<< DPeakToCompton << endl;
132  txtfile << "Peak/Compton (res):\t"<< PeakToCompton << " +- "<< DPeakToCompton << endl;
133  }
134  else
135  {
136 
137  cout << "SumCompton = 0 -> no PeakToCompton possible!"<<endl;
138  txtfile << "SumCompton = 0 -> no PeakToCompton possible!"<<endl;
139  }
140  Double_t AllEntries = hGammaSpectrum->GetEntries();
141  cout << "Entries: " << AllEntries << endl;
142  txtfile << "Entries: " << AllEntries << endl;
143 
144 
145  //fitting the histograms
146 
147  TF1 *GausBG = new TF1("GausBG","gausn"/*(0)+pol1(3)" */,(PeakX-100)/ChannelResolution,(PeakX+100)/ChannelResolution);
148  GausBG->SetParName(0,"Ampl");
149  GausBG->SetParName(1,"x0");
150  GausBG->SetParName(2,"sigma");
151  GausBG->SetParameter(0,PeakY/10000);
152  GausBG->SetParameter(1,PeakX/Double_t(ChannelResolution));
153  GausBG->SetParameter(2,100./Double_t(ChannelResolution));
154  GausBG->SetParLimits(0,0.0000001, 1000000);
155  GausBG->SetParLimits(1,(PeakX-10)/Double_t(ChannelResolution),(PeakX+10)/Double_t(ChannelResolution));
156  GausBG->SetParLimits(2,1./Double_t(ChannelResolution),100./Double_t(ChannelResolution));
157  /*GausBG->SetParameter(0,1);
158  GausBG->SetParameter(1,1);
159  GausBG->SetParameter(2,0);
160  GausBG->FixParameter(2,0);
161  GausBG->FixParameter(1,0);
162  GausBG->FixParameter(0,0);
163  */
164  GausBG->SetNpx(100000);
165  GausBG->SetLineColor(kRed);
166  hGammaSpectrum->Fit(GausBG,"R");
167 
168  txtfile << "Fitparameter Spectrum (model: gausn):" << endl;
169  for (Int_t i = 0; i < 3; i++)
170  {
171  txtfile << GausBG->GetParName(i) << ":\t" << GausBG->GetParameter(i) << endl;
172  }
173  Double_t FWHM = GausBG->GetParameter(2)*2.3548200;
174  txtfile << "FWHM[keV]:\t" << FWHM*1000000 << endl;
175  cout << "FWHM[keV]:\t" << FWHM*1000000 << endl;
176 
177  /*TF1 *Poisson = new TF1("Poisson","gausn",0,15);//,2); //maybe gausn? maybe no fit? //gausn --> no ,2);
178  //Poisson->SetParameter(0,4);
179  //Poisson->SetParameter(1,SumPeak);
180  //Poisson->SetParLimits(0,0,12);
181  Poisson->SetParName(0,"Ampl");//Poisson->SetParName(0,"lambda");
182  Poisson->SetParName(1,"x0");//Poisson->SetParName(1,"Ampl");
183  Poisson->SetParName(2,"sigma");
184  //Poisson->SetParameter(1,1000);
185  Poisson->SetNpx(100000);
186  Poisson->SetLineColor(kRed);
187  //hNoHits->Fit(Poisson);*/
188 
189  //txtfile << "Fitparameter Number of Hits (model: Ampl.*Poisson):" << endl;
190  //for (Int_t i = 0; i < 3; i++) // if gausn fit --> < 2 -> < 3 , if no fit, comment it out!!!
191  //{
192  // txtfile << Poisson->GetParName(i) << ":\t" << Poisson->GetParameter(i) << endl;
193  //}
194 
195  //some make-up for the histograms
196  Double_t FullEnergyPeakEff = double(int(GausBG->GetParameter(0)*ChannelResolution/(nEvents*2)*100*1000))/1000;
197  Double_t FullEnergyPeakEffError = double(int(GausBG->GetParError(0)*ChannelResolution/(nEvents*2)*100*1000))/1000;
198 
199  cout << "Full-Energy-Peak-Eff. [%]: " << FullEnergyPeakEff << endl;; // *2 because only 2Pi of solid angle is simulated
200  txtfile << "Full-Energy-Peak-Eff. [%]:\t" << FullEnergyPeakEff << endl; // *2 because only 2Pi of solid angle is simulated
201  // writing to files and closing
202  cout << "Error ofFull-Energy-Peak-Eff. [%]: " << FullEnergyPeakEffError<< endl; // *2 because only 2Pi of solid angle is simulated
203  txtfile << "Error of Full-Energy-Peak-Eff. [%]:\t" << FullEnergyPeakEffError << endl; // *2 because only 2Pi of solid angle is simulated
204  // writing to files and closing
205  TFile *Output = new TFile("OutputCombinedAna.root","RECREATE");
206  hGammaSpectrum->Write();
207  hNoOfHits->Write();
208 
209  txtfile.close();
210  Output->Close();
211  File->Close();
212 
213  TString CollectFileName = Anapath +"/CollectFEPF.txt";
214  ofstream CollectFile(CollectFileName, ios::out|ios::app); //collect name and Full Energy Eff of all Sims here
215  CollectFile << Geometry << "\t" << Sektar << "\t"<< Energy << "\t"<< nEvents << "\t" << FullEnergyPeakEff << "\t" << FullEnergyPeakEffError << "\t" <<FWHM*1000000 << "\t"<< PeakToCompton<< "\t"<<InputName << endl;
216 
217  CollectFile.close();
218 
219  cout << "macro finished successfully" << endl;
220 
221  return 0;
222 }
i
Int_t i
Definition: run_full.C:25
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
s
TLorentzVector s
Definition: Pnd2DStar.C:50
FileName
TString FileName
Definition: modifyRootFscGeometryFile.C:269
invCompton
Double_t invCompton(Double_t Energy, Double_t Eprime)
Definition: AnalyseCombinedGammaFiles.C:22
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
nEvents
Int_t nEvents
Definition: hit_dirc.C:11
out
TFile * out
Definition: reco_muo.C:20
Compton
Double_t Compton(Double_t Energy, Double_t Th)
Definition: AnalyseCombinedGammaFiles.C:18
can3
TCanvas * can3
Definition: GammaSpectraAnalysis.C:89
Double_t Compton ( Double_t  Energy,
Double_t  Th 
)

Definition at line 18 of file AnalyseCombinedGammaFiles.C.

References CAMath::Cos().

Referenced by AnalyseCombinedGammaFiles(), GammaSpectraAnalysis_CableTest(), and GammaSpectraAnalysis_NoH().

19 {
20 return Energy-Energy/(1+Energy/0.000511*(1-TMath::Cos(Th)));
21 }
CAMath::Cos
static T Cos(const T &x)
Definition: PndCAMath.h:43
Double_t invCompton ( Double_t  Energy,
Double_t  Eprime 
)

Definition at line 22 of file AnalyseCombinedGammaFiles.C.

Referenced by AnalyseCombinedGammaFiles(), GammaSpectraAnalysis_CableTest(), and GammaSpectraAnalysis_NoH().

23 {
24 return TMath::ACos((Energy*Energy-Energy*Eprime-Eprime*0.000511)/(Energy*(Energy-Eprime)));
25 }
Double_t PeakFunc ( Double_t *  x,
Double_t *  par 
)

Definition at line 26 of file AnalyseCombinedGammaFiles.C.

References Pi, and CAMath::Sqrt().

27 {
28  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]);
29 }
CAMath::Sqrt
static T Sqrt(const T &x)
Definition: PndCAMath.h:37
par
Double_t par[3]
Definition: anaLmdCluster.C:146
x
Double_t x
Definition: createRootFscGeometryFile.C:159
Pi
Double_t Pi
Definition: PndGeoHypGeCrystal.cxx:10
Double_t PoissonFunc ( Double_t *  x,
Double_t *  par 
)

Definition at line 30 of file AnalyseCombinedGammaFiles.C.

31 {
32  return par[1]*TMath::Poisson(x[0],par[0]);
33 }
par
Double_t par[3]
Definition: anaLmdCluster.C:146
x
Double_t x
Definition: createRootFscGeometryFile.C:159
Generated on Wed Apr 3 2019 10:03:10 for FairRoot/PandaRoot by   doxygen 1.8.5