#include <PndHypGeGammaAna.h>

Inheritance diagram for PndHypGeGammaAna:

Public Member Functions
	PndHypGeGammaAna ()

	PndHypGeGammaAna (TString TxtFileName_Ext, Double_t GammaEnergy_Ext, Int_t nEvents=0, Int_t nPeaks_Ext=1, Int_t PeakToLook_Ext=1)

	~PndHypGeGammaAna ()

virtual InitStatus	Init ()

virtual void	Exec (Option_t *opt)

void	Finish ()

void	SetParContainers ()

void	SetTarget (Bool_t TargetAvailable=1)

void	SetOmegaQuadrupolMode (Bool_t useOmegaQuadrupolMode_ext=1)

void	SetPeakWidtchStrechFactor (Double_t PeakWidtchStrechFactor_ext=1)

Protected Member Functions
Double_t	Compton (Double_t E, Double_t Th)

Double_t	invCompton (Double_t E, Double_t Eprime)

Double_t	SmearHit (Double_t Energy)

Double_t	CalculatePeakWidth (Double_t Energy)

void	FitSpectrum ()

void	FittingOutput ()

void	FindPeaks ()

void	GetPeakContent ()

void	ActivateComptonCalculations (Bool_t isActivated=1)

void	CalculateCompton ()

void	WriteHistogramsToFile ()

void	HistogramCosmetics (TH1D *histo, TString XTitle, Double_t XTitleSize, Double_t XTitleOffset, Double_t XLabelSize, TString YTitle, Double_t YTitleSize, Double_t YTitleOffset, Double_t YLabelSize)

void	HistogramCosmetics2D (TH2D *histo, TString XTitle, Double_t XTitleSize, Double_t XTitleOffset, Double_t XLabelSize, TString YTitle, Double_t YTitleSize, Double_t YTitleOffset, Double_t YLabelSize)

bool	InTargetAbsorbedGammaWouldHitGermanium ()

	ClassDef (PndHypGeGammaAna, 1)

Protected Attributes
TString	TxtFileName

ofstream	TxtFile

Bool_t	useOmegaQuadrupolMode

Double_t	ResConstPar

Double_t	ResGradientPar

Double_t	PeakWidtchStrechFactor

Double_t	PeakWidth

Double_t	bufferEnergy

TH1D *	hGamEnergy

TH1D *	hGamEnergyTargetBackground

TH1D *	hNumberOfHits

TH1D *	hThetaCheck

TH1D *	hThetaCheckPrimary

TH1D *	hAbsorption

TH2D *	h2AbsorptionDistanceAngle

TClonesArray *	fMcTr

TClonesArray *	fMc

TClonesArray *	fGe

TClonesArray *	fGeAl

TClonesArray *	fSilicon

TClonesArray *	fTargetOther

PndHypGePoint *	fHitGe

PndHypPoint *	fHitSi

PndHypPoint *	fHitTargetOther

PndMCTrack *	fMCGam

PndMCTrack *	fMCGamPrimary

Long_t	EvtCount

Long_t	NumberOfEvents

Double_t	GammaEnergy

Double_t	iHistoUpperLimit

TString	fName

TRandom *	Rng

FairRootManager *	ioman

std::set< int >	SetOfCrystalHit

TVector3	fStartVertex

TString	fVertexVolumeName

Int_t	nPeaks

Int_t	PeakToLook

Double_t	PeakX

Double_t	PeakY

Double_t	DPeakY

Double_t	SumPeak

Double_t	DSumPeak

Bool_t	DoComptonCalculations

Double_t	lowComptonAngle

Double_t	highComptonAngle

Double_t	SumCompton

Double_t	DSumCompton

Double_t	PeakToCompton

Double_t	DPeakToCompton

TF1 *	GausBG

Float_t *	xpeaks

Float_t *	ypeaks

Bool_t	fTargetAvailable

Detailed Description

Definition at line 27 of file PndHypGeGammaAna.h.

Constructor & Destructor Documentation

PndHypGeGammaAna::PndHypGeGammaAna ( )

Definition at line 23 of file PndHypGeGammaAna.cxx.

References DoComptonCalculations, and fTargetAvailable.

         : FairTask("Panda HypGeGammaAna  Task"), xpeaks(0), ypeaks(0), fHitGe(0), fMCGam(0), fHitSi(0), fHitTargetOther(0)
 {
         DoComptonCalculations = 0;
         fTargetAvailable=0;
 
 }

PndHypGeGammaAna::PndHypGeGammaAna	(	TString	TxtFileName_Ext,
		Double_t	GammaEnergy_Ext,
		Int_t	nEvents = `0`,
		Int_t	nPeaks_Ext = `1`,
		Int_t	PeakToLook_Ext = `1`
	)

Definition at line 30 of file PndHypGeGammaAna.cxx.

References DoComptonCalculations, fTargetAvailable, GammaEnergy, nEvents, nPeaks, NumberOfEvents, PeakToLook, PeakWidtchStrechFactor, TxtFileName, and useOmegaQuadrupolMode.

         : FairTask("Panda HypGeGammaAna  Task"), xpeaks(0), ypeaks(0), fHitGe(0), fMCGam(0), fHitSi(0), fHitTargetOther(0)
 {
         TxtFileName=TxtFileName_Ext;
         GammaEnergy=GammaEnergy_Ext;
         NumberOfEvents= nEvents;
         nPeaks= nPeaks_Ext;
         PeakToLook = PeakToLook_Ext;
         DoComptonCalculations = 0;
         fTargetAvailable=0;
         useOmegaQuadrupolMode=0;
         PeakWidtchStrechFactor = 1;
 }

PndHypGeGammaAna::~PndHypGeGammaAna ( )

Definition at line 44 of file PndHypGeGammaAna.cxx.

References fGe, fGeAl, fHitGe, fHitSi, fHitTargetOther, fMc, fMCGam, fMcTr, fSilicon, fTargetOther, hAbsorption, hGamEnergy, hNumberOfHits, Rng, xpeaks, and ypeaks.

 {
         //delete fStorage;
         delete hGamEnergy;
         delete hNumberOfHits;
         delete hAbsorption;
         delete Rng;
         delete xpeaks;
         delete ypeaks;
         delete fHitGe;
         delete fMCGam;
         delete fHitSi;
         delete fHitTargetOther;
 
         if (fMcTr)
         {
                 fMcTr->Delete();
                 delete fMcTr;
         }
         if (fMc)
         {
                 fMc->Delete();
                 delete fMc;
         }
         if (fGe)
         {
                 fGe->Delete();
                 delete fGe;
         }
         if (fGeAl)
         {
                 fGeAl->Delete();
                 delete fGeAl;
         }
         if (fSilicon)
         {
                 fSilicon->Delete();
                 delete fSilicon;
         }
         if (fTargetOther)
         {
                 fTargetOther->Delete();
                 delete fTargetOther;
         }
 }

Member Function Documentation

void PndHypGeGammaAna::ActivateComptonCalculations ( Bool_t isActivated = 1 )

protected

Definition at line 566 of file PndHypGeGammaAna.cxx.

References DoComptonCalculations.

 {
         DoComptonCalculations= isActivated;
 }

void PndHypGeGammaAna::CalculateCompton ( )

protected

Definition at line 571 of file PndHypGeGammaAna.cxx.

References Compton(), DPeakToCompton, DPeakY, DSumCompton, fVerbose, GammaEnergy, hGamEnergy, highComptonAngle, i, lowComptonAngle, PeakToCompton, PeakY, sqrt(), SumCompton, and TxtFile.

 {
 // get number of compton events
 
         for (Int_t i =Compton(GammaEnergy,lowComptonAngle)*1000*100000; i <= Compton(GammaEnergy,highComptonAngle)*1000*100000; i++)
   {
                 SumCompton += hGamEnergy->GetBinContent(i);
         }
         DSumCompton = sqrt(SumCompton);
 
         if (fVerbose >= 2)
         {
                 cout <<"Hyp Ge Gamma Ana:\tCompton from bins [" << Compton(GammaEnergy,lowComptonAngle)*1000*100000 << " to "<< Compton(GammaEnergy,highComptonAngle)*1000*100000 << "]" << endl;
                 cout <<"Hyp Ge Gamma Ana:\tNumber of compton events: "<< SumCompton<< " +- " << DSumCompton << endl;
         }
         TxtFile <<"Hyp Ge Gamma Ana:\tCompton range [bins] \t[" << Compton(GammaEnergy,lowComptonAngle)*1000*100000 << " to "<< Compton(GammaEnergy,highComptonAngle)*1000*100000 << "]" << endl;
         TxtFile <<"Hyp Ge Gamma Ana:\tNumber of compton events: \t"<< SumCompton<< " +- " << DSumCompton << endl;
 
         // calculate Peak/Compton ratio
 
         if (SumCompton != 0)
         {
                 PeakToCompton = PeakY/SumCompton*(Compton(GammaEnergy,highComptonAngle)-Compton(GammaEnergy,lowComptonAngle))*1000*100000;
                 DPeakToCompton = sqrt(pow(DPeakY/SumCompton,2)+ pow(PeakY*DSumCompton/SumCompton/SumCompton,2))*(Compton(GammaEnergy,highComptonAngle)-Compton(GammaEnergy,lowComptonAngle))*1000*100000;
                 if (fVerbose >= 2)
                         cout << "Hyp Ge Gamma Ana:\tPeak/Compton : "<< PeakToCompton << " ± "<< DPeakToCompton << endl;
                 TxtFile << "Hyp Ge Gamma Ana:\tPeak/Compton : \t"<< PeakToCompton << " ± "<< DPeakToCompton << endl;
         }
         else
         {
                 if (fVerbose >= 2)
                         cout << "Hyp Ge Gamma Ana:\tSumCompton = 0 -> no PeakToCompton possible!"<<endl;
                 TxtFile << "Hyp Ge Gamma Ana:\tSumCompton = 0 -> no PeakToCompton possible!"<<endl;
         }
 }

Double_t PndHypGeGammaAna::CalculatePeakWidth ( Double_t Energy )

protected

Definition at line 440 of file PndHypGeGammaAna.cxx.

References PeakWidtchStrechFactor, ResConstPar, ResGradientPar, and sqrt().

Referenced by Init(), and SmearHit().

 {
         return PeakWidtchStrechFactor *sqrt(ResConstPar + Energy * ResGradientPar);
 }

PndHypGeGammaAna::ClassDef	(	PndHypGeGammaAna	,
		1
	)

protected

Double_t PndHypGeGammaAna::Compton	(	Double_t	E,
		Double_t	Th
	)

protected

Definition at line 424 of file PndHypGeGammaAna.cxx.

References CAMath::Cos().

Referenced by CalculateCompton().

 {
 return E-E/(1+E/0.000511*(1-TMath::Cos(Th)));
 }

void PndHypGeGammaAna::Exec ( Option_t * opt )

virtual

Definition at line 208 of file PndHypGeGammaAna.cxx.

References CAMath::Abs(), Double_t, EvtCount, fGe, fHitGe, fHitTargetOther, fMCGam, fMCGamPrimary, fMcTr, fSilicon, fStartVertex, fTargetAvailable, fTargetOther, fVerbose, fVertexVolumeName, GammaEnergy, PndHypGePoint::GetDetectorID(), PndHypPoint::GetDetName(), PndHypGePoint::GetEnergyLoss(), PndMCTrack::GetMomentum(), PndMCTrack::GetMotherID(), PndMCTrack::GetStartVertex(), PndHypGePoint::GetTrackID(), PndHypPoint::GetXin(), PndHypPoint::GetZin(), gGeoManager, h2AbsorptionDistanceAngle, hAbsorption, hGamEnergy, hGamEnergyTargetBackground, hNumberOfHits, hThetaCheck, hThetaCheckPrimary, i, NumberOfEvents, PeakWidth, Pi, SetOfCrystalHit, SmearHit(), and TString.

 {
         Double_t EnergyPerCrystal[48] = {0};
         Double_t CollectedEnergy = 0;
         SetOfCrystalHit.clear();
         bool isFirstInteraction = 0;
         if (fTargetAvailable)
         {
                 if( fSilicon->GetEntriesFast() || fTargetOther->GetEntriesFast())
                 {
                         bool isEndl = 0;
                         fMCGamPrimary = (PndMCTrack*)fMcTr->At(0);
                         TVector3 PrimaryVertex = fMCGamPrimary->GetStartVertex();
                         TVector3 PrimaryMomentum = fMCGamPrimary->GetMomentum();
                         Double_t PrimaryDistance;
                         Double_t PrimaryTheta;
 //                      if (fSilicon->GetEntriesFast())
 //                      {
 //                              fHitSi = (PndHypPoint*)fSilicon->At(0);
 //                              fMCGam = (PndMCTrack*)fMcTr->At(fHitSi->GetTrackID());
 //                              if (fHitSi->GetTrackID() == 1 && fMCGam->GetMotherID() <1)
 //                              {
 //                                      isFirstInteraction = 1;
 //                                      //cout << "TrPIDSi\t\t" << fMCGam->GetPdgCode()<<"\t\t";
 //                                      isEndl = 1;
 //                                      fStartVertex = fMCGam->GetStartVertex();
 //                                      PrimaryDistance = (fStartVertex - PrimaryVertex).Mag();
 //                                      PrimaryTheta = PrimaryMomentum.Theta()*180/TMath::Pi();
 //                                      fVertexVolumeName = fHitSi->GetDetName();//=gGeoManager->FindNode(fStartVertex.X(),fStartVertex.Y(),fStartVertex.Z())->GetName();
 //                                      if(PrimaryTheta >110 && PrimaryTheta <165) // range of germaniums
 //                                                      hAbsorption->Fill(1);
 //                              }
 //                      }
                         if (fTargetOther->GetEntriesFast())
                         {
                                 fHitTargetOther = (PndHypPoint*)fTargetOther->At(0);
                                 fMCGam = (PndMCTrack*)fMcTr->At(fHitTargetOther->GetTrackID());
                                 if (fHitTargetOther->GetTrackID() == 1 && fMCGam->GetMotherID() < 1 )
                                 {
                                         //isFirstInteraction = 1;
                                         //cout << "TrPIDOt\t\t" << fMCGam->GetPdgCode()<<"\t\t";
                                         isEndl = 1;
                                         fStartVertex = fMCGam->GetStartVertex();
                                         TVector3 InVector(fHitTargetOther->GetXin(),fHitTargetOther->GetZin(),fHitTargetOther->GetZin());
                                         TVector3 InDiffVector =fStartVertex - InVector;
                                         TVector3 DiffVector =fStartVertex - PrimaryVertex;
                                         PrimaryDistance = DiffVector.Mag();
                                         PrimaryTheta = DiffVector.Theta()*180/TMath::Pi();
                                         fVertexVolumeName = fHitTargetOther->GetDetName();//=gGeoManager->FindNode(fStartVertex.X(),fStartVertex.Y(),fStartVertex.Z())->GetName();
                                         TString StartVolumeName =gGeoManager->FindNode(fStartVertex.X(),fStartVertex.Y(),fStartVertex.Z())->GetName();
                                         if(PrimaryTheta >110 && PrimaryTheta <165 && (StartVolumeName.Contains("Sensor") ||StartVolumeName.Contains("Absorber") ||StartVolumeName.Contains("HoldingFrame") ||StartVolumeName.Contains("Window") ||StartVolumeName.Contains("Frame"))) // range of germaniums
                                         {
                                                 isFirstInteraction = 1;
                                                 if (PrimaryDistance >10)
                                                 {
                                                         cout << PrimaryDistance << endl;
                                                         cout << PrimaryTheta << endl;
 
                                                         cout << fVertexVolumeName.Data() << endl;
                                                         cout << StartVolumeName.Data() << endl;
                                                         fStartVertex.Print();
                                                         PrimaryVertex.Print();
                                                         DiffVector.Print();
                                                         //InDiffVector.Print();
                                                         cout << fHitTargetOther->GetLength() << endl;
                                                 }
                                                 if (fVertexVolumeName.Contains("Sensor"))
                                                 {
                                                         hAbsorption->Fill(1);
                                                 }
                                                 else
                                                 {
                                                         if (fVertexVolumeName.Contains("Absorber"))
                                                         {
                                                                 hAbsorption->Fill(2);
                                                                 //cout << fVertexVolumeName.Data()<< endl;
                                                         }
                                                         else
                                                         {
                                                                 if (fVertexVolumeName.Contains("HoldingFrame"))
                                                                         hAbsorption->Fill(3);
                                                                 else
                                                                 {
                                                                         if (fVertexVolumeName.Contains("Window"))
                                                                                 hAbsorption->Fill(4);
                                                                         else
                                                                         {
                                                                                 if (fVertexVolumeName.Contains("Frame"))
                                                                                         hAbsorption->Fill(5);
                                                                                 else
                                                                                 {
                                                                                         hAbsorption->Fill(6);
                                                                                         cout << fVertexVolumeName.Data()<< endl;
                                                                                         fStartVertex.Print();
                                                                                 }
                                                                         }
                                                                 }
                                                         }
                                                 }
                                         }
                                 }
                         }       // end of if (fTargetOther->GetEntriesFast())
                         if(isFirstInteraction)
                         {
                                 if(PrimaryTheta >110 && PrimaryTheta <165)
                                         h2AbsorptionDistanceAngle->Fill(PrimaryTheta,PrimaryDistance);
                         }
                         //if (isEndl)
                                 //cout << endl;
                 }       // end of if( fSilicon->GetEntriesFast() || fTargetOther->GetEntriesFast())
         }               // end of if (fTargetAvailable)
         // loop over hits and save energy loss per crystal and crystals hit
 
         // theta angle of gammas hiting the germanium
         fMCGam = (PndMCTrack*)fMcTr->At(0);
         TVector3 GammaMomentum = fMCGam->GetMomentum();
         //bool ishThetaCheckPrimaryFilled=0;
         if (fGe->GetEntriesFast())
         {
                 hThetaCheck->Fill(GammaMomentum.Theta()*180/TMath::Pi());
                 if (fTargetAvailable)
                 {
                         if (!(fSilicon->GetEntriesFast() || fTargetOther->GetEntriesFast()))
                                 hThetaCheckPrimary->Fill(GammaMomentum.Theta()*180/TMath::Pi());
                 }
                 else
                 {
                         hThetaCheckPrimary->Fill(GammaMomentum.Theta()*180/TMath::Pi());
                         //ishThetaCheckPrimaryFilled=1;
                 }
         }
         //if (!ishThetaCheckPrimaryFilled && fGe->GetEntriesFast()&& fTargetAvailable && !(fSilicon->GetEntriesFast() || fTargetOther->GetEntriesFast()))
         //{
                 //hThetaCheck->Fill(GammaMomentum.Theta()*180/TMath::Pi());
                 //if (GammaMomentum.Theta()*180/TMath::Pi() > 170)
                 //{
                         //cout << GammaMomentum.Theta()*180/TMath::Pi() << endl;
                         //cout << fMCGam->GetMotherID() << endl;
                 //}
         //}
         // gamma energy part
         for (Int_t i=0; i < fGe->GetEntriesFast(); i++)
         {
                 fHitGe=(PndHypGePoint*)fGe->At(i);
                 fMCGam = (PndMCTrack*)fMcTr->At(fHitGe->GetTrackID());
                 if (fHitGe->GetEnergyLoss() !=0)
                 {
                         EnergyPerCrystal[fHitGe->GetDetectorID() % 100] += fHitGe->GetEnergyLoss()*1000;                        // in MeV
                         SetOfCrystalHit.insert(fHitGe->GetDetectorID() % 100);
                 }
         }
 
         // add up the whole collected energy by all crystals, smear every crystal seperatly
         for (std::set<int>::iterator it=SetOfCrystalHit.begin(); it!=SetOfCrystalHit.end(); ++it)
         {
                 CollectedEnergy += SmearHit(EnergyPerCrystal[*it]);
         }
         // fill histogram if there is any collected Energy
         if (CollectedEnergy > 0)
         {
                 hGamEnergy->Fill(CollectedEnergy);
                 if (fTargetAvailable)
                 {
                         if (fSilicon->GetEntriesFast() || fTargetOther->GetEntriesFast())
                                 hGamEnergyTargetBackground->Fill(CollectedEnergy);
                 }
                 if (fVerbose >= 2 && CollectedEnergy != 0)
                         cout << "Hyp Ge Gamma Ana:\t" << "Energy " << CollectedEnergy << endl;
         }
         // fill number of hits histo when full gamma energy is collected, NoH = size of the set
         if (TMath::Abs(CollectedEnergy-GammaEnergy) < 10*PeakWidth)                     //Peak = 10 * PeakWidth of 1 crystal --> takes multiple hits into account (faster than fitting and than doing it again, error is small)
         {                       
                 hNumberOfHits->Fill(SetOfCrystalHit.size());            //Fill # of Hits -diagramm
                 if (fVerbose >= 2)
                         cout << "Hyp Ge Gamma Ana:\t" << "NoHFill: " << SetOfCrystalHit.size() << endl;
         }
         // show progress every 10 % of events
         if (!((EvtCount++*10) % NumberOfEvents))
         {
                 if (EvtCount == 1)
                         cout << "Hyp Ge Gamma Ana:\t"<<"0  % done!" << endl;
                 else
                         cout << "Hyp Ge Gamma Ana:\t" << 100 * Double_t(EvtCount-1)/NumberOfEvents << " % done!" << endl;
         }
         if (fVerbose == 3)
                 cout << "Hyp Ge Gamma Ana:\t" << "Event Number"<< EvtCount << endl;
 }

void PndHypGeGammaAna::FindPeaks ( )

protected

Definition at line 523 of file PndHypGeGammaAna.cxx.

References DPeakY, fVerbose, hGamEnergy, i, nPeaks, PeakToLook, PeakX, PeakY, s, sqrt(), xpeaks, and ypeaks.

 {
         TSpectrum *s = new TSpectrum(nPeaks);
         s->Search(hGamEnergy,0.01,"goffnodraw",0.01);
 
         if (fVerbose >= 2)
                         cout << "Hyp Ge Gamma Ana:\tSpectrum search finished: " << hGamEnergy->GetEntries() <<endl;
 
         xpeaks = s->GetPositionX();
         ypeaks = s->GetPositionY();
 
         if (fVerbose == 3)
         {
                 for (Int_t i = 0; i < nPeaks; i++)
                 {
                         cout <<"Hyp Ge Gamma Ana:\tPeak @ ("<< xpeaks[i]<<"," << ypeaks[i] << ")" << endl;
                 }
         }
         PeakX = xpeaks[PeakToLook-1]*100000;
         PeakY = ypeaks[PeakToLook-1];
         DPeakY = sqrt(PeakY);
 
         delete s;
 }

void PndHypGeGammaAna::Finish ( )

Definition at line 396 of file PndHypGeGammaAna.cxx.

References fVerbose, hGamEnergy, TxtFile, and WriteHistogramsToFile().

 {
         //-------------Analysis of the spectrum-------------------------------------
         cout << "Hyp Ge Gamma Ana:\tEvent loop finished" <<endl;
         //print number of entries in spectrum
         if (fVerbose >= 2)
                 cout << "Hyp Ge Gamma Ana:\tEntries in spectrum: " << hGamEnergy->GetEntries() <<endl;
         TxtFile << "Hyp Ge Gamma Ana:\tEntries in spectrum: \t" << hGamEnergy->GetEntries() <<endl;
 
 
         //FindPeaks();
 
         //if (DoComptonCalculations)
         //      CalculateCompton();
 
         
         //FitSpectrum();
         //GetPeakContent();
         //FittingOutput();
                 //--------------make up for the histograms----------------------------------
 
         WriteHistogramsToFile();
 
         TxtFile.close();
         cout << "Hyp Ge Gamma Ana:\tAnalysis finished succesfully" << endl;
 }

void PndHypGeGammaAna::FitSpectrum ( )

protected

Definition at line 469 of file PndHypGeGammaAna.cxx.

References GausBG, hGamEnergy, PeakX, and PeakY.

 {
         //fitting the spectrum                                          modell is not good!!!! M.Steinen
         GausBG = new TF1("GausBG","gausn",(PeakX-1000)/100000,(PeakX+1000)/100000);
         GausBG->SetParName(0,"Ampl");
         GausBG->SetParName(1,"x0");
         GausBG->SetParName(2,"sigma");
         GausBG->SetParameter(0,PeakY/10000);
         GausBG->SetParameter(1,PeakX/100000.);
         GausBG->SetParameter(2,100./100000.);
         GausBG->SetParLimits(0,0.000001, 1000000);
         GausBG->SetParLimits(1,(PeakX-100)/100000.,(PeakX+100)/100000.);
         GausBG->SetParLimits(2,1./100000.,1000./100000.);
         /*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);
         hGamEnergy->Fit(GausBG,"R");
 }

void PndHypGeGammaAna::FittingOutput ( )

protected

Definition at line 494 of file PndHypGeGammaAna.cxx.

References Double_t, fVerbose, GausBG, i, NumberOfEvents, and TxtFile.

 {
 cout << "Hyp Ge Gamma Ana:\tCHi²red of the fit \t"<<GausBG->GetChisquare() / GausBG->GetNDF() << endl;
         cout << "Hyp Ge Gamma Ana:\tCHi²red of the fit \t"<<GausBG->GetChisquare() <<"  "<< GausBG->GetNDF() << endl;
         TxtFile << "Hyp Ge Gamma Ana:\tCHi²red of the fit \t"<<GausBG->GetChisquare() / GausBG->GetNDF() << endl;
         TxtFile << "Hyp Ge Gamma Ana:\tFitparameter Spectrum (model: gausn):" << endl;
         for (Int_t i = 0; i < 3; i++)
         {
                 TxtFile <<"Hyp Ge Gamma Ana:\t" << GausBG->GetParName(i) << ":\t" << GausBG->GetParameter(i) << endl;
         }
         if (fVerbose >= 2)
         {
                 cout << "Hyp Ge Gamma Ana:\tFitparameter Spectrum (model: gausn):" << endl;
                 for (Int_t i = 0; i < 3; i++)
                 {
                         cout <<"Hyp Ge Gamma Ana:\t" << GausBG->GetParName(i) << ":\t" << GausBG->GetParameter(i) << endl;
                 }
         }
 
         Double_t FWHM = GausBG->GetParameter(2)*2.3548200;
         if (fVerbose >= 2)
                 cout << "Hyp Ge Gamma Ana:\tFWHM[keV]: " <<  FWHM*1000 << endl;
         TxtFile << "Hyp Ge Gamma Ana:\tFWHM[keV]:\t" <<  FWHM*1000 << endl;
 
         if (fVerbose >= 2)
                 cout << "Hyp Ge Gamma Ana:\tFull-Energy-Peak-Eff. [%]: " << double(int(GausBG->GetParameter(0)*100000/(NumberOfEvents*2)*100*1000))/1000 << " ± " << double(int(GausBG->GetParError(0)*100000/(NumberOfEvents*2)*100*1000))/1000 << endl; // *2 because only 2Pi of solid angle is simulated
         TxtFile << "Hyp Ge Gamma Ana:\tFull-Energy-Peak-Eff. [%]: \t" << double(int(GausBG->GetParameter(0)*100000/(NumberOfEvents*2)*100*1000))/1000 << " ± " << double(int(GausBG->GetParError(0)*100000/(NumberOfEvents*2)*100*1000))/1000 << endl; // *2 because only 2Pi of solid angle is simulated
 }

void PndHypGeGammaAna::GetPeakContent ( )

protected

Definition at line 548 of file PndHypGeGammaAna.cxx.

References DPeakY, DSumPeak, fVerbose, hGamEnergy, i, PeakWidth, PeakX, PeakY, sqrt(), SumPeak, and TxtFile.

 {
         //get content of peak
         SumPeak = 0;
         for(Int_t i = PeakX-PeakWidth*10*1000000*100; i <= PeakX+PeakWidth*10*1000000*100; i++)         //Peak = 10 * PeakWidth of 1 crystal --> takes multiple hits into account (faster than fitting and than doing it again, error is small)
                 SumPeak += hGamEnergy->GetBinContent(i);
         DSumPeak = sqrt(SumPeak);
         if (fVerbose >= 2)
         {
                 cout <<"Hyp Ge Gamma Ana:\tPeak at bin "<< PeakX << endl;
                 cout <<"Hyp Ge Gamma Ana:\tPeak bin height"<< PeakY << " ± " << DPeakY<< endl;
                 cout <<"Hyp Ge Gamma Ana:\tSumPeak "<<SumPeak  << " ± " << DSumPeak << endl;
         }
         TxtFile <<"Hyp Ge Gamma Ana:\tPeak at bin \t"<< PeakX << endl;
         TxtFile <<"Hyp Ge Gamma Ana:\tPeak bin height\t"<< PeakY << " ± " << DPeakY<< endl;
         TxtFile <<"Hyp Ge Gamma Ana:\tSumPeak \t"<<SumPeak << " ± " << DSumPeak << endl;
 }

void PndHypGeGammaAna::HistogramCosmetics	(	TH1D *	histo,
		TString	XTitle,
		Double_t	XTitleSize,
		Double_t	XTitleOffset,
		Double_t	XLabelSize,
		TString	YTitle,
		Double_t	YTitleSize,
		Double_t	YTitleOffset,
		Double_t	YLabelSize
	)

protected

Definition at line 445 of file PndHypGeGammaAna.cxx.

Referenced by WriteHistogramsToFile().

 {
         histo->SetXTitle(XTitle.Data());
         histo->GetXaxis()->SetTitleSize(XTitleSize);
         histo->GetXaxis()->SetTitleOffset(XTitleOffset);
         histo->GetXaxis()->SetLabelSize(XLabelSize);
 
         histo->SetYTitle(YTitle.Data());
         histo->GetYaxis()->SetTitleSize(YTitleSize);
         histo->GetYaxis()->SetTitleOffset(YTitleOffset);
         histo->GetYaxis()->SetLabelSize(YLabelSize);
 }

void PndHypGeGammaAna::HistogramCosmetics2D	(	TH2D *	histo,
		TString	XTitle,
		Double_t	XTitleSize,
		Double_t	XTitleOffset,
		Double_t	XLabelSize,
		TString	YTitle,
		Double_t	YTitleSize,
		Double_t	YTitleOffset,
		Double_t	YLabelSize
	)

protected

Definition at line 457 of file PndHypGeGammaAna.cxx.

Referenced by WriteHistogramsToFile().

 {
         histo->SetXTitle(XTitle.Data());
         histo->GetXaxis()->SetTitleSize(XTitleSize);
         histo->GetXaxis()->SetTitleOffset(XTitleOffset);
         histo->GetXaxis()->SetLabelSize(XLabelSize);
 
         histo->SetYTitle(YTitle.Data());
         histo->GetYaxis()->SetTitleSize(YTitleSize);
         histo->GetYaxis()->SetTitleOffset(YTitleOffset);
         histo->GetYaxis()->SetLabelSize(YLabelSize);
 }

InitStatus PndHypGeGammaAna::Init ( )

virtual

Definition at line 91 of file PndHypGeGammaAna.cxx.

References CalculatePeakWidth(), DoComptonCalculations, EvtCount, fGe, fGeAl, fMcTr, fName, fSilicon, fTargetOther, fVerbose, GammaEnergy, gGeoManager, h2AbsorptionDistanceAngle, hAbsorption, hGamEnergy, hGamEnergyTargetBackground, highComptonAngle, hNumberOfHits, hThetaCheck, hThetaCheckPrimary, iHistoUpperLimit, invCompton(), ioman, lowComptonAngle, NumberOfEvents, PeakWidth, ResConstPar, ResGradientPar, Rng, TxtFile, TxtFileName, and useOmegaQuadrupolMode.

 {
         cout << "Hyp Ge Gamma Ana:\t" << "Starting Init of HypGe Gamma Ana" << endl;
         
         //----------get FairRootManager-------------------------------------
         ioman = FairRootManager::Instance();
   if ( ! ioman ) {
     cout << "-E- PndEmcHitProducer::Init: "
          << "RootManager not instantiated!" << endl;
     return kFATAL;
   }
    
         if (fVerbose == 3)
                 cout << "Hyp Ge Gamma Ana:\t" << "Getting FairRootManager finished" << endl;
                 
         //----------get input arrays---------------------------------------------
         fMcTr = (TClonesArray*) ioman->GetObject("MCTrack");
         fGe = (TClonesArray*) ioman->GetObject("HypGePoint");
         fGeAl = (TClonesArray*) ioman->GetObject("HypGeAlPoint");
         fSilicon = (TClonesArray*) ioman->GetObject("HypPoint");
         fTargetOther = (TClonesArray*) ioman->GetObject("HypSTMatBudPoint");
         
 
 
         if (fVerbose ==3)
                 cout << "Hyp Ge Gamma Ana:\t" << "Getting input arrays finished" << endl;
                 
         //-----------extracting file name--------------------------------------
         fName = ioman->GetInFile()->GetName();
         if (fVerbose == 3)
         {
                 cout << "Hyp Ge Gamma Ana:\t" << "Full input file name with path: " << fName << endl;
                 cout << "Hyp Ge Gamma Ana:\t" << "Chopping first " << fName.Last('/')+1 << " chars"<< endl;
         }
         fName.Remove(0,fName.Last('/')+1);
         if (fVerbose == 3)
                 cout << "Hyp Ge Gamma Ana:\t" << "File name: " <<fName<< endl;
 
         if (fVerbose == 3)
                 cout << "Hyp Ge Gamma Ana:\t" << "File name extraction finished" << endl;
                 
         //----------create histograms------------------------------------------
         if (useOmegaQuadrupolMode)
         {
                 iHistoUpperLimit = 0.6;
         }
         else
         {
                 iHistoUpperLimit = GammaEnergy*1.2;
         }
         hGamEnergy = new TH1D("GamEnergy","Total gamma energy deposit, " + fName, 10000*iHistoUpperLimit,0.0,iHistoUpperLimit);
         hGamEnergyTargetBackground = new TH1D("hGamEnergyTargetBackground","gamma background due to reactions in target, " + fName, 10000*iHistoUpperLimit,0.0,iHistoUpperLimit);
         hNumberOfHits= new TH1D("NumberOfHits", "Number of Hits, " + fName, 16,0,15);
         hThetaCheck= new TH1D("hThetaCheck", "Theta Check, " + fName, 90,90,180);
         hThetaCheckPrimary= new TH1D("hThetaCheckPrimary", "Theta Check only primaries, " + fName, 90,90,180);
         
         hAbsorption = new TH1D("hAbsorption", "Part of Gamma interaction in target system," + fName, 10,1,11);
                 hAbsorption ->GetXaxis()->SetBinLabel(1,"Silicon");                                                     // gamma interaction in si detectors
                 hAbsorption ->GetXaxis()->SetBinLabel(2,"Absorber");                                            // gamma interaction in absorber
                 hAbsorption ->GetXaxis()->SetBinLabel(3,"HoldingFrame");                                                                // gamma interaction in holding frame
                 hAbsorption ->GetXaxis()->SetBinLabel(4,"Window foil");                                 // gamma interaction in window foil
                 hAbsorption ->GetXaxis()->SetBinLabel(5,"beam pipe");                                                           // gamma interaction in beam pipe
                 hAbsorption ->GetXaxis()->SetBinLabel(6,"Other");                                                               // gamma interaction in other parts
         h2AbsorptionDistanceAngle = new TH2D("h2AbsorptionDistanceAngle","Correlation of absorption distance and angle," + fName,45,90,180, 100,0,50);
 
         if (fVerbose == 3)
                 cout << "Hyp Ge Gamma Ana:\t" << "Histogram creation finished" << endl;
         
         //----------set detector resolution---------------------------------------------
                         //TO DO: more realistic resolution (smearing) modell    M.Steinen
         ResConstPar =1.40872e-07;
         ResGradientPar = 4.16193e-07 ;
         PeakWidth = CalculatePeakWidth(GammaEnergy);
         //----------initialize EvtCount----------------------------------------
         EvtCount = 0;
         
         //----------initialize Compton angles-----------------------
         if (DoComptonCalculations)
         {
                 lowComptonAngle = invCompton(0.001332,0.001040);
                 highComptonAngle = invCompton(0.001332,0.001096);
         }
 
         //----------Get number of events from FairRootManager-----------------------
         
         if (NumberOfEvents== 0)
                 NumberOfEvents= ioman->GetInChain()->GetEntries();
         cout << "Hyp Ge Gamma Ana:\tNumber of analyzed events: "<< NumberOfEvents << endl;
                 
         //----------initialize random generator------------------------------------
         Rng = new TRandom(time(NULL));
         cout << "Hyp Ge Gamma Ana:\t" << "Init of HypGe Gamma Ana finished successfully" << endl;
         
         //----------create storage in root file------------------------------------
         //fStorage = new PndHypGeGammaAnaStorage("AnaData","Analysis Data");
         //      fStorage-> SetAllEntries(NumberOfEvents);
                                                 //fStorage-> SetInFile((ioman->GetInFile()->GetName()).Data());                 // this line doen't work!!!!! why?!?!?!?!?!??!?!
         //----------open txt outputfile--------------------------------------------
 
         TxtFile.open(TxtFileName);
                 TxtFile << "Hyp Ge Gamma Ana:\tFile to analyze: \t"<< ioman->GetInFile()->GetName() << endl;
                 TxtFile << "Hyp Ge Gamma Ana:\tNumber of analyzed events: \t"<< NumberOfEvents << endl;
 
 
         std::cout << "-I- gGeoManager = "<<gGeoManager << std::endl;
         return kSUCCESS;
 }

bool PndHypGeGammaAna::InTargetAbsorbedGammaWouldHitGermanium ( )

protected

Definition at line 648 of file PndHypGeGammaAna.cxx.

 {
         //to be implemented
         // this function checks if a gamma that was absorbed in the target system would have hit the solid angle of the germanium, only than it's important for the material budget
         // to do this,  the start vertex and the momentum vector of the primary gamma have to be projected and checked of this line crosses a germanium detector
         return 0;
 }

Double_t PndHypGeGammaAna::invCompton	(	Double_t	E,
		Double_t	Eprime
	)

protected

Definition at line 428 of file PndHypGeGammaAna.cxx.

Referenced by Init().

 {
 return TMath::ACos((E*E-E*Eprime-Eprime*0.000511)/(E*(E-Eprime)));
 }

void PndHypGeGammaAna::SetOmegaQuadrupolMode ( Bool_t useOmegaQuadrupolMode_ext = 1 )

Definition at line 660 of file PndHypGeGammaAna.cxx.

References useOmegaQuadrupolMode.

Referenced by GammaSpectraAnalysis_NoH_Task().

 {
         useOmegaQuadrupolMode=useOmegaQuadrupolMode_ext;
 }

void PndHypGeGammaAna::SetParContainers ( )

Definition at line 199 of file PndHypGeGammaAna.cxx.

References rtdb.

 {
   // Get Base Container
   FairRun* ana = FairRun::Instance();
   FairRuntimeDb* rtdb=ana->GetRuntimeDb();
   std::cout << "-I- rtdb = "<<rtdb << std::endl;
 
 }

void PndHypGeGammaAna::SetPeakWidtchStrechFactor ( Double_t PeakWidtchStrechFactor_ext = 1 )

Definition at line 665 of file PndHypGeGammaAna.cxx.

References PeakWidtchStrechFactor.

Referenced by GammaSpectraAnalysis_NoH_Task().

 {
         PeakWidtchStrechFactor=PeakWidtchStrechFactor_ext;
 }

void PndHypGeGammaAna::SetTarget ( Bool_t TargetAvailable = 1 )

Definition at line 656 of file PndHypGeGammaAna.cxx.

References fTargetAvailable.

Referenced by GammaSpectraAnalysis_NoH_Task().

 {
         fTargetAvailable=TargetAvailable;
 }

Double_t PndHypGeGammaAna::SmearHit ( Double_t Energy )

protected

Definition at line 433 of file PndHypGeGammaAna.cxx.

References CalculatePeakWidth(), Double_t, and Rng.

Referenced by Exec().

 {
         Double_t Resolution =  CalculatePeakWidth(Energy);
         Double_t SmearedEnergy = Energy +Rng->Gaus(0,Resolution);
         return SmearedEnergy;
 }

void PndHypGeGammaAna::WriteHistogramsToFile ( )

protected

Definition at line 607 of file PndHypGeGammaAna.cxx.

References h2AbsorptionDistanceAngle, hAbsorption, hGamEnergy, hGamEnergyTargetBackground, HistogramCosmetics(), HistogramCosmetics2D(), hNumberOfHits, hThetaCheck, and hThetaCheckPrimary.

Referenced by Finish().

 {
 //-------------writing and drawing the histograms-------------------------
 
 HistogramCosmetics(hGamEnergy, "#gamma Energy [MeV]", 0.06, 0.7, 0.045, "Counts", 0.06, 0.8, 0.045);
 HistogramCosmetics(hGamEnergyTargetBackground, "#gamma Energy [MeV]", 0.06, 0.7, 0.045, "Counts", 0.06, 0.8, 0.045);
 HistogramCosmetics(hNumberOfHits, "Number of Crystals hit", 0.06, 0.7, 0.045, "Counts", 0.06, 0.8, 0.045);
 HistogramCosmetics(hThetaCheck, "#Theta [#circ]", 0.06, 0.7, 0.045, "Counts", 0.06, 0.8, 0.045);
 HistogramCosmetics(hThetaCheckPrimary, "#Theta [#circ]", 0.06, 0.7, 0.045, "Counts", 0.06, 0.8, 0.045);
 HistogramCosmetics(hAbsorption, "part of absorption", 0.06, 0.7, 0.045, "Counts", 0.06, 0.8, 0.045);
 HistogramCosmetics2D(h2AbsorptionDistanceAngle, "#Theta [#circ]", 0.06, 0.7, 0.045, "Distance [cm]", 0.06, 0.8, 0.045);
 
 hGamEnergy->Write();
 hGamEnergyTargetBackground->Write();
 hNumberOfHits->Write();
 hThetaCheck->Write();
 hThetaCheckPrimary->Write();
 hAbsorption->Write();
 h2AbsorptionDistanceAngle->Write();
 TCanvas* CanGamEnergy = new TCanvas("CanGamEnergy", "Gamma Energy", 1000,600);
 hGamEnergy->Draw();
 
 TCanvas* CanGamEnergyTarget = new TCanvas("CanGamEnergyTarget", "Gamma Energy target", 1000,600);
 hGamEnergyTargetBackground->Draw();
 
 TCanvas* CanNoH = new TCanvas("CanNoH", "Number of Hits", 1000,600);
 hNumberOfHits->Draw();
 
 TCanvas* CanTH = new TCanvas("CanTH", "TH", 1000,600);
 hThetaCheck->Draw();
 
 TCanvas* CanTHPrimary = new TCanvas("CanTHPrimary", "TH prim", 1000,600);
 hThetaCheckPrimary->Draw();
 
 TCanvas* CanAbsorption = new TCanvas("CanAbsorption", "part of absorption", 1000,600);
 hAbsorption->Draw();
 
 TCanvas* CanAbsorption2D = new TCanvas("CanAbsorption2D", "correlation of absorption distance and theta", 1000,600);
 h2AbsorptionDistanceAngle->Draw("colz");
 }