/******************************************************

Analysis Task created by M.Steinen steinen@kph.uni-mainz.de
Analysis of Gamma Simulation with hypGe detectors
*******************************************************/ 


#include "PndHypGeGammaAna.h"
#include <iostream>
#include <TMath.h>
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include <time.h>
#include <TCanvas.h>
#include <TSpectrum.h>

//#include "RhoBase/TFactory.h"


using namespace std;


PndHypGeGammaAna::PndHypGeGammaAna()
	: 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, Int_t nPeaks_Ext,Int_t PeakToLook_Ext) 
	: 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()
{
	//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;
	}
}


InitStatus PndHypGeGammaAna::Init()
{
	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;
}

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

}

void PndHypGeGammaAna::Exec(Option_t*)
{
	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::Finish()
{
	//-------------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;
}


Double_t PndHypGeGammaAna::Compton(Double_t E, Double_t Th)
{
return E-E/(1+E/0.000511*(1-TMath::Cos(Th)));
}
Double_t PndHypGeGammaAna::invCompton(Double_t E, Double_t Eprime)
{
return TMath::ACos((E*E-E*Eprime-Eprime*0.000511)/(E*(E-Eprime)));
}

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

Double_t PndHypGeGammaAna::CalculatePeakWidth(Double_t Energy)
{
	return PeakWidtchStrechFactor *sqrt(ResConstPar + Energy * ResGradientPar);
}

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)
{
	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)
{
	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::FitSpectrum()
{
	//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()
{
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::FindPeaks()
{
	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::GetPeakContent()
{
	//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::ActivateComptonCalculations(Bool_t isActivated)
{
	DoComptonCalculations= isActivated;
}

void PndHypGeGammaAna::CalculateCompton()
{
// 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;
	}
}

void PndHypGeGammaAna::WriteHistogramsToFile()
{
//-------------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");
}

bool PndHypGeGammaAna::InTargetAbsorbedGammaWouldHitGermanium()
{
	//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;
}

void PndHypGeGammaAna::SetTarget(Bool_t  TargetAvailable)
{
	fTargetAvailable=TargetAvailable;
}
void PndHypGeGammaAna::SetOmegaQuadrupolMode(Bool_t useOmegaQuadrupolMode_ext )
{
	useOmegaQuadrupolMode=useOmegaQuadrupolMode_ext;
}

void PndHypGeGammaAna::SetPeakWidtchStrechFactor(Double_t PeakWidtchStrechFactor_ext)
{
	PeakWidtchStrechFactor=PeakWidtchStrechFactor_ext;
}

ClassImp(PndHypGeGammaAna)