// root macro to analyze the simulation output
//void convertMCPoints()
{

  // -----  Load libraries   ------------------------------------------------
//``gSystem->Load("fstream.h");
   gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");

  
   gSystem->Load("libHypGe");
  
  // -----   Timer   --------------------------------------------------------
  TStopwatch timer;
  timer.Start();
  // ------------------------------------------------------------------------
#include <vector>




	// the sim file you want to analyse
	string Filename = "../sim_pseudodata/sim_hypGe_BgUrqmd2_Marcell_Co60_3er_highStat"; //without File Type ending!!!
	string CompleteFilename = Filename+".root";
  TFile* g = new TFile(CompleteFilename.c_str());
	
//Output Files
 	TString outfile= "../Data_Marcell/"+Filename+".root";
  TFile* fi = new TFile(outfile,"RECREATE");
	TString txtfileName ="../Data_Marcell/"+Filename+".txt";
	ofstream txtfile;
	txtfile.open(txtfileName,ios::trunc);
	txtfile << "File read:" << CompleteFilename << endl;
   //photons from hyp electromag. decay
  TTree *b=(TTree *) g->Get("cbmsim") ;
  TClonesArray* hit_bar=new TClonesArray("PndHypGePoint");
  b->SetBranchAddress("HypGePoint",&hit_bar);//Branch names
  TClonesArray* mc_bar=new TClonesArray("PndMCTrack");
  b->SetBranchAddress("MCTrack",&mc_bar);//Branch names


  //**** Energy deposited from particle background
      
  
  //****photons from hyp elect. decay
	double upperborder = 0.009;
  string Name = "total gam energy deposit, " + Filename;
  TH1D* gamTde = new TH1D("gamTde",Name.c_str(),10000000*upperborder,0.0000,upperborder);
 
  
	bool verbose = false;
	Int_t MotherId,Motherpdg;
	
	
	
	TVector3 vecs,pos;
	int mcpdg = -1,ev;
	Double_t mult,En,Eng,Enth;
	
	//vector<int> event;
	int count;
	cout<<b->GetEntriesFast()<<endl;
	
Double_t Resolution = 2;	//keV
	for (Int_t k=0; k<b->GetEntriesFast(); k++)
	  { 
			cout << k << endl;
	    Eng=0.;
	    b->GetEntry(k);
	    //if(verbose) cout<<"Event No "<<j<<endl;
	    for (Int_t i=0; i<hit_bar->GetEntriesFast(); i++)
	      { 
					cout << hit_bar->GetEntriesFast()<<endl;
					PndHypGePoint *hitgam=(PndHypGePoint*)hit_bar->At(i);
		
					PndMCTrack *mcgam = (PndMCTrack*)mc_bar->At(hitgam->GetTrackID());
		
					Double_t rande;
					//rande= gRandom->Gaus(hitgam->GetEnergyLoss(),0.000003);
					//if (sci->Getpdgcode())
		
					Eng +=gRandom->Gaus(0,0.000001*Resolution/2.3548);
	
					//if (sci)cout<<sci->GetEnergyLoss()<<endl;
					//Eng +=gRandom->Gaus(0,0.000001*Resolution);
					//randy= gRandom->Gaus(0,1);
					
		
	
		
					Eng =Eng + (hitgam->GetEnergyLoss());
					//Eng =Eng + (rande);
		
	      }//end for i (points in event)
	    //count =0;
			//cout <<Eng<<endl;
			
	    if(Eng>0)gamTde->Fill(Eng);
	   
	    
	    
	  }// end for j (events)
  

  
	TCanvas* can3 = new TCanvas("can3","germanium detector",0,0,1000,1000);
 
  gamTde->Draw();
Double_t PeakX = 40000, PeakY = gamTde->GetBinContent(13320), DPeakY = sqrt(PeakY); 
Double_t SumPeak = 0;
for(int i = PeakX-Resolution*10*3/2.3548; i <= PeakX+Resolution*10*3/2.3548; i++)
	SumPeak += gamTde->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 SumCompton=0, DSumCompton = 0; 
  for (int i =9623; i <= 11176; i++)
	SumCompton += gamTde->GetBinContent(i);
DSumCompton = sqrt(SumCompton);

cout <<"N. of Compton Events: "<< SumCompton<< " +- " << DSumCompton << endl;
txtfile <<"N. of Compton Events: "<< SumCompton<< " +- " << DSumCompton << endl;
//cout << "Peak/Compton: "<< gamTde->GetMaximum()/SumCompton << " +- "<< sqrt(pow(DPeakY/SumCompton,2)+ pow(PeakY*DSumCompton/SumCompton/SumCompton,2)) << endl;
  
txtfile << "Peak/Compton (res): "<< SumPeak/SumCompton << " +- "<< sqrt(pow(DSumPeak/SumCompton,2)+ pow(SumPeak*DSumCompton/SumCompton/SumCompton,2)) << endl;   
Double_t AllEntries = gamTde->GetEntries();
cout << "Entries: " << AllEntries << endl;
txtfile << "Entries: " << AllEntries << endl;
gamTde->Write();
fi->Close();
txtfile.close();
	// -----   Finish   -------------------------------------------------------
	timer.Stop();
	Double_t rtime = timer.RealTime();
	Double_t ctime = timer.CpuTime();
	cout << endl << endl;
	cout << "Macro finished succesfully." << endl;
	cout << "Output file is "    << outfile << endl;
	cout << "Parameter file is " << txtfileName << endl;
	cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	cout << endl;
	// ------------------------------------------------------------------------
	
}