//----------------------------------------------------------------------
// File and Version Information:
//      $Id: //
// Description:
//      Class PndEmcMcTruthWriter. Module to take the ADC waveforms and produces digi.
// 
//	 Software developed for the BaBar Detector at the SLAC B-Factory.
// Adapted for the PANDA experiment at GSI		
//		
// Author List:
//      Phil Strother                  Original Author
// Dima Melnichuk - adaption for PANDA		
// Copyright Information:
//      Copyright (C) 1996             Imperial College
//
//----------------------------------------------------------------------

#include "PndEmcMcTruthWriter.h"
#include "PndEmcCluster.h"
#include "PndEmcCorrection.h"
#include "FairEventHeader.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "PndEmcBump.h"
#include "TClonesArray.h"
#include "TStopwatch.h"
#include <vector>		
#include <iostream>
#include "PndMCTrack.h"
#include "PndEmcGeoPar.h"
#include "PndEmcStructure.h"

//#include <map>
#include <functional>
//#include "TLorentzVector.h"

using std::cout;
using std::endl;
using std::fstream;

/*struct key 
{
	key(Int_t hit, Double_t e) : iHit(hit), fEnergy(e){}
	Int_t iHit;
	Double_t fEnergy;
};
bool operator< (const key& lhs, const key& rhs) {
	return lhs.fEnergy < rhs.fEnergy;
}
bool operator> (const key& lhs, const key& rhs) {
	return lhs.fEnergy > rhs.fEnergy;
}*/

PndEmcMcTruthWriter::PndEmcMcTruthWriter(Int_t verbose, Bool_t storeHist):
	/*fRecoHitArray(0),*/ fStoreHist(storeHist), fResult(0), tResult(0)
{
}
/*ostream& operator<<(ostream& os, const TVector3& mom) {
	os<<"("<<mom.Px()<<", "<<mom.Py()<<", "<<mom.Pz()<<")";
	os<<std::endl;
	return os;
}*/

//--------------
// Destructor --
//--------------

PndEmcMcTruthWriter::~PndEmcMcTruthWriter()
{
}


InitStatus PndEmcMcTruthWriter::Init()
{
	// Get RootManager
	FairRootManager* ioman = FairRootManager::Instance();
	if ( ! ioman )
	{
		cout << "-E- PndEmcMcTruthWriter::Init: "
			<< "RootManager not instantiated!" << endl;
		return kFATAL;
	}

	// Get input array
/*	fRecoHitArray = (TClonesArray*) ioman->GetObject("EmcBump");
	if ( ! fRecoHitArray ) {
		cout << "-W- PndEmcMcTruthWriter::Init: "
			<< "No PndEmcBump array!" << endl;
		return kERROR;
	}

	fCluCorrArray = (TClonesArray*) ioman->GetObject("EmcBumpCorr");
	if ( ! fRecoHitArray ) {
		cout << "-W- PndEmcMcTruthWriter::Init: "
			<< "No PndEmcBumpCorr array!" << endl;
		return kERROR;
	}*/

	fMCTrackArray = (TClonesArray*) ioman->GetObject("MCTrack");
	if ( ! fMCTrackArray ) {
		cout << "-W- PndEmcMakeCluster::Init: "
			<< "No MCTrack array! Needed for MC Truth" << endl;
		//    return kERROR;
	}

	fResult = new TFile("gammaeta_truth.root", "RECREATE");
	tResult = new TTree("Res1", "result");
	//tResult->Branch("E1", &fEnergy1, "E1/D");
	//tResult->Branch("E1C", &fEnergy1C, "E1C/D");
	//tResult->Branch("E1CC", &fEnergy1CC, "E1CC/D");
	//tResult->Branch("p1", &p4_1[0], "p1[4]/D");
	tResult->Branch("Mcp1", &Mcp4_1[0], "Mcp1[4]/D");
	//tResult->Branch("v1", &v_1[0], "v1[3]/D");
	tResult->Branch("Theta1", &fTheta1, "Theta1/D");
	tResult->Branch("Phi1", &fPhi1, "Phi1/D");
	tResult->Branch("X1", &fX1, "X1/D");
	tResult->Branch("Y1", &fY1, "Y1/D");
	tResult->Branch("Z1", &fZ1, "Z1/D");
	//tResult->Branch("Z201", &fZ201, "Z201/D");
	//tResult->Branch("Z531", &fZ532, "Z531/D");
	//tResult->Branch("Lat1", &fLat1, "Lat1/D");
	//tResult->Branch("Mod1", &fModule1, "Mod1/I");
	//tResult->Branch("McE1", &fMCE1, "McE1/D");
	//tResult->Branch("McPx1", &fMCPx1, "McPx1/D");
	//tResult->Branch("McPy1", &fMCPy1, "McPy1/D");
	//tResult->Branch("McPz1", &fMCPz1, "McPz1/D");
	//tResult->Branch("McVx1", &fMCVx1, "McVx1/D");
	//tResult->Branch("McVy1", &fMCVy1, "McVy1/D");
	//tResult->Branch("McVz1", &fMCVz1, "McVz1/D");
	//tResult->Branch("McT1", &fMCT1, "McT1/D");

	//tResult->Branch("E2", &fEnergy2, "E2/D");
	//tResult->Branch("E2C", &fEnergy2C, "E2C/D");
	//tResult->Branch("E2CC", &fEnergy2CC, "E2CC/D");
	//tResult->Branch("p2", &p4_2[0], "p2[4]/D");
	tResult->Branch("Mcp2", &Mcp4_2[0], "Mcp2[4]/D");
	//tResult->Branch("v2", &v_2[0], "v2[3]/D");
	tResult->Branch("Theta2", &fTheta2, "Theta2/D");
	tResult->Branch("Phi2", &fPhi2, "Phi2/D");
	tResult->Branch("X2", &fX2, "X2/D");
	tResult->Branch("Y2", &fY2, "Y2/D");
	tResult->Branch("Z2", &fZ2, "Z2/D");
	tResult->Branch("Det1", &fDet1, "Det1/I");
	tResult->Branch("Det2", &fDet2, "Det2/I");
	tResult->Branch("Det3", &fDet3, "Det3/I");
	//tResult->Branch("Z202", &fZ202, "Z202/D");
	//tResult->Branch("Z532", &fZ532, "Z532/D");
	//tResult->Branch("Lat2", &fLat2, "Lat2/D");
	//tResult->Branch("Mod2", &fModule2, "Mod2/I");
	//tResult->Branch("McE2", &fMCE2, "McE2/D");
	//tResult->Branch("McPx2", &fMCPx2, "McPx2/D");
	//tResult->Branch("McPy2", &fMCPy2, "McPy2/D");
	//tResult->Branch("McPz2", &fMCPz2, "McPz2/D");
	//tResult->Branch("McVx2", &fMCVx2, "McVx2/D");
	//tResult->Branch("McVy2", &fMCVy2, "McVy2/D");
	//tResult->Branch("McVz2", &fMCVz2, "McVz2/D");
	//tResult->Branch("McT2", &fMCT2, "McT2/D");

	//tResult->Branch("E3", &fEnergy3, "E3/D");
	//tResult->Branch("E3C", &fEnergy3C, "E3C/D");
	//tResult->Branch("E3CC", &fEnergy3CC, "E3CC/D");
	//tResult->Branch("p3", &p4_3[0], "p3[4]/D");
	tResult->Branch("Mcp3", &Mcp4_3[0], "Mcp3[4]/D");
	//tResult->Branch("v3", &v_3[0], "v3[3]/D");
	tResult->Branch("Theta3", &fTheta3, "Theta3/D");
	tResult->Branch("Phi3", &fPhi3, "Phi3/D");
	tResult->Branch("X3", &fX3, "X3/D");
	tResult->Branch("Y3", &fY3, "Y3/D");
	tResult->Branch("Z3", &fZ3, "Z3/D");
	tResult->Branch("isHit1", &fHitEmc1, "isHit1/I");
	tResult->Branch("isHit2", &fHitEmc2, "isHit2/I");
	tResult->Branch("isHit3", &fHitEmc3, "isHit3/I");
	//tResult->Branch("Z203", &fZ203, "Z203/D");
	//tResult->Branch("Z533", &fZ533, "Z533/D");
	//tResult->Branch("Lat3", &fLat3, "Lat3/D");
	//tResult->Branch("Mod3", &fModule3, "Mod3/I");
	//tResult->Branch("m12", &fm12, "m12/D");
	//tResult->Branch("m23", &fm23, "m23/D");
	//tResult->Branch("m31", &fm31, "m31/D");
	//tResult->Branch("m123", &fm123, "m123/D");
	//tResult->Branch("MCm12", &fMCm12, "MCm12/D");
	//tResult->Branch("MCm23", &fMCm23, "MCm23/D");
	//tResult->Branch("MCm31", &fMCm31, "MCm31/D");
	//tResult->Branch("MCm123", &fMCm123, "MCm123/D");
	//tResult->Branch("McE3", &fMCE3, "McE3/D");
	//tResult->Branch("McPx3", &fMCPx3, "McPx3/D");
	//tResult->Branch("McPy3", &fMCPy3, "McPy3/D");
	//tResult->Branch("McPz3", &fMCPz3, "McPz3/D");
	//tResult->Branch("McVx3", &fMCVx3, "McVx3/D");
	//tResult->Branch("McVy3", &fMCVy3, "McVy3/D");
	//tResult->Branch("McVz3", &fMCVz3, "McVz3/D");
	//tResult->Branch("McT3", &fMCT3, "McT3/D");

	//tResult->Branch("ExtE", &fOtherEnergy, "ExtE/D");
	//tResult->Branch("ExtEC", &fOtherEnergyC, "ExtEC/D");
	//tResult->Branch("NDigi", &fNumDigi, "NDigi/I");
	//tResult->Branch("NBump", &fNumBump, "NBump/I");
	//tResult->Branch("pJpsi", &fJpsi[0], "pJpsi[4]/D");
	//tResult->Branch("pMcJpsi", &fMCJpsi[0], "pMcJpsi[4]/D");


	//fGeoPar->InitEmcMapper();
	//fMapper=PndEmcMapper::Instance();
	//fEmcStr=PndEmcStructure::Instance();

	cout << "-I- PndEmcMcTruthWriter: Intialization successfull" << endl;

	return kSUCCESS;
}

void PndEmcMcTruthWriter::Exec(Option_t* opt)
{
	TStopwatch timer;
	if (fVerbose>2){
		timer.Start();
	}

	//fOtherEnergy = 0.;
	//fOtherEnergyC = 0.;
	Double_t EventTime = FairRootManager::Instance()->GetEventTime();

	//Int_t nHits = fRecoHitArray->GetEntriesFast();

	static Int_t evtNo = 0;
	cout<<"**************************************"<<endl;
	//Int_t evtNo = FairRunAna::Instance()->GetEventHeader()->GetMCEntryNumber();
	cout<<"Event No. #"<<(evtNo)<<", EvtTime #"<<EventTime<<std::endl;
	//cout<<"PndEmcMcTruthWriter: "<<nHits<<" photon found."<<endl;
	cout<<"**************************************"<<endl;
	//std::vector<key> goodHits;
	//std::vector<Int_t> theMostEnergy(0,3);
	//if(nHits >= 3){
	//	for (Int_t iHit = 0; iHit <nHits; iHit++) {
	//		PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(iHit);
	//		Double_t Energy = theHit->energy();
	//		if(Energy < 0.10) continue;//100 MeV 
	//		goodHits.push_back(key(iHit, Energy));
	//	}

	//	bool photon1_found(false);
	//	bool photon2_found(false);
	//	bool photon3_found(false);
	//	if(goodHits.size() >=3)
	//	{
	//		sort(goodHits.begin(), goodHits.end(), std::greater<key>());
	//		Double_t angle(999.);
	//		Int_t HitIndex1;
	//		for(Int_t i=0;i<goodHits.size();++i){
	//			PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(goodHits[i].iHit);
	//			TVector3 hit_mom = theHit->where();
	//			PndMCTrack* pt1 =((PndMCTrack*)fMCTrackArray->At(0));
	//			Double_t angle1 = hit_mom.Angle(pt1->GetMomentum());
	//			//Double_t angle1 = fabs(theHit->GetEnergyCorrected() - pt1->Get4Momentum().E());
	//			if( angle1 < angle){
	//				angle = angle1;
	//				HitIndex1 = goodHits[i].iHit;
	//			}
	//		}
	//		std::cout<<"Find Photon 1 #"<<angle<<std::endl;
	//		{
	//			PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(HitIndex1);
	//			fEnergy1 = theHit->energy();
	//			fEnergy1C = theHit->GetEnergyCorrected()/1.0165;
	//			PndEmcCorrection* theCorr = (PndEmcCorrection*)fCluCorrArray->At(HitIndex1);
	//			fEnergy1CC = theCorr->EnergyCorrPhoton();
	//			//TVector3 pos = theHit->position();
	//			v_1 = theHit->where();
	//			fTheta1 = v_1.Theta();//theCorr->ThetaCorrPhoton();//v_1.Theta();
	//			fPhi1 = v_1.Phi();
	//			p4_1 = TLorentzVector(fEnergy1C*TMath::Sin(fTheta1)*TMath::Cos(fPhi1)
	//					, fEnergy1C*TMath::Sin(fTheta1)*TMath::Sin(fPhi1)
	//					, fEnergy1C*TMath::Cos(fTheta1), fEnergy1C);
	//			//fX1 = theHit->x();
	//			//fY1 = theHit->y();
	//			//fZ1 = theHit->z();
	//			fZ201 = theHit->Z20();
	//			fZ531 = theHit->Z53();
	//			fLat1 = theHit->LatMom();
	//			fModule1 = theHit->GetModule();
	//		}
	//		angle = 999.;
	//		Int_t HitIndex2;
	//		for(Int_t i=0;i<goodHits.size();++i){
	//			PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(goodHits[i].iHit);
	//			TVector3 hit_mom = theHit->where();
	//			PndMCTrack* pt2 = ((PndMCTrack*)fMCTrackArray->At(1));
	//			Double_t angle2 = hit_mom.Angle(pt2->GetMomentum());
	//			//Double_t angle2 = fabs(theHit->GetEnergyCorrected() - pt2->Get4Momentum().E());
	//			if( (goodHits[i].iHit != HitIndex1) && angle2 < angle){
	//				angle = angle2;
	//				HitIndex2 = goodHits[i].iHit;
	//			}
	//		}
	//		std::cout<<"Find Photon 2 #"<<angle<<std::endl;
	//		{
	//			PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(HitIndex2);
	//			fEnergy2 = theHit->energy();
	//			fEnergy2C = theHit->GetEnergyCorrected()/1.0165;
	//			PndEmcCorrection* theCorr = (PndEmcCorrection*)fCluCorrArray->At(HitIndex2);
	//			fEnergy2CC = theCorr->EnergyCorrPhoton();
	//			v_2 = theHit->where();
	//			fTheta2 = v_2.Theta();//theCorr->ThetaCorrPhoton();//v_2.Theta();
	//			fPhi2 = v_2.Phi();
	//			p4_2 = TLorentzVector(fEnergy2C*TMath::Sin(fTheta2)*TMath::Cos(fPhi2)
	//					, fEnergy2C*TMath::Sin(fTheta2)*TMath::Sin(fPhi2)
	//					, fEnergy2C*TMath::Cos(fTheta2), fEnergy2C);
	//			//fX2 = theHit->x();
	//			//fY2 = theHit->y();
	//			//fZ2 = theHit->z();
	//			fZ202 = theHit->Z20();
	//			fZ532 = theHit->Z53();
	//			fLat2 = theHit->LatMom();
	//			fModule2 = theHit->GetModule();
	//		}
	//		angle = 999.;
	//		Int_t HitIndex3;
	//		for(Int_t i=0;i<goodHits.size();++i){
	//			PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(goodHits[i].iHit);
	//			TVector3 hit_mom = theHit->where();
	//			PndMCTrack* pt3 = ((PndMCTrack*)fMCTrackArray->At(3));
	//			Double_t angle3 = hit_mom.Angle(pt3->GetMomentum());
	//			//Double_t angle3 = fabs(theHit->GetEnergyCorrected() - pt3->Get4Momentum().E());
	//			if( (goodHits[i].iHit != HitIndex1 && goodHits[i].iHit != HitIndex2) && angle3 < angle){
	//				angle = angle3;
	//				HitIndex3 = goodHits[i].iHit;
	//			}
	//		}
	//		std::cout<<"Find Photon 3 #"<<angle<<std::endl;
	//		{
	//			PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(HitIndex3);
	//			fEnergy3 = theHit->energy();
	//			fEnergy3C = theHit->GetEnergyCorrected()/1.0165;
	//			PndEmcCorrection* theCorr = (PndEmcCorrection*)fCluCorrArray->At(HitIndex3);
	//			fEnergy3CC = theCorr->EnergyCorrPhoton();
	//			v_3 = theHit->where();
	//			fTheta3 = v_3.Theta();//theCorr->ThetaCorrPhoton();//v_3.Theta();
	//			fPhi3 = v_3.Phi();
	//			p4_3 = TLorentzVector(fEnergy3C*TMath::Sin(fTheta3)*TMath::Cos(fPhi3)
	//					, fEnergy3C*TMath::Sin(fTheta3)*TMath::Sin(fPhi3)
	//					, fEnergy3C*TMath::Cos(fTheta3), fEnergy3C);
	//			//fX3 = theHit->x();
	//			//fY3 = theHit->y();
	//			//fZ3 = theHit->z();
	//			fZ203 = theHit->Z20();
	//			fZ533 = theHit->Z53();
	//			fLat3 = theHit->LatMom();
	//			fModule3 = theHit->GetModule();
	//		}
	//		for(Int_t i=0;i<goodHits.size();++i){
	//			if(goodHits[i].iHit != HitIndex1 && goodHits[i].iHit != HitIndex2 && goodHits[i].iHit != HitIndex3){
	//				PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(goodHits[i].iHit);
	//				fOtherEnergy += theHit->energy();
	//				fOtherEnergyC += theHit->GetEnergyCorrected();
	//			}
	//		}
	//		std::cout<<"The most energetic photon #"<<goodHits[0].iHit<<", "<<goodHits[1].iHit<<", "<<goodHits[2].iHit<<std::endl;
	//		std::cout<<"Angle match #"<<HitIndex1 <<", "<<HitIndex2 <<", "<<HitIndex3<<std::endl;


	//		fNumBump = goodHits.size();
	//		//fNumDigi = theHit->NumberOfDigis();
	//		PndMCTrack* pt1 =((PndMCTrack*)fMCTrackArray->At(0));
	//		Mcp4_1 = TLorentzVector(pt1->GetMomentum(), pt1->GetE());
	//		//fMCE1 =  pt1->GetE();
	//		//fMCPx1 = pt1->GetMomentum().Px();
	//		//fMCPy1 = pt1->GetMomentum().Py();
	//		//fMCPz1 = pt1->GetMomentum().Pz();
	//		//fMCVx1 = pt1->GetStartVertex().X();
	//		//fMCVy1 = pt1->GetStartVertex().Y();
	//		//fMCVz1 = pt1->GetStartVertex().Z();
	//		//fMCT1 =  pt1->GetStartTime();
	//		//cout<<"p1 #"<<fMCPx1<<", "<<fMCPy1<<", "<<fMCPz1<<", "<<fMCE1<<std::endl;
	//		Mcp4_1.Print();

	//		PndMCTrack* pt2 =((PndMCTrack*)fMCTrackArray->At(1));
	//		Mcp4_2 = TLorentzVector(pt2->GetMomentum(), pt2->GetE());
	//		//fMCE2 =  pt2->GetE();
	//		//fMCPx2 = pt2->GetMomentum().Px();
	//		//fMCPy2 = pt2->GetMomentum().Py();
	//		//fMCPz2 = pt2->GetMomentum().Pz();
	//		//fMCVx2 = pt2->GetStartVertex().X();
	//		//fMCVy2 = pt2->GetStartVertex().Y();
	//		//fMCVz2 = pt2->GetStartVertex().Z();
	//		//fMCT2 =  pt2->GetStartTime();
	//		//cout<<"p2 #"<<fMCPx2<<", "<<fMCPy2<<", "<<fMCPz2<<", "<<fMCE2<<std::endl;
	//		Mcp4_2.Print();

	//		PndMCTrack* pt3 =((PndMCTrack*)fMCTrackArray->At(2));
	//		Mcp4_3 = TLorentzVector(pt3->GetMomentum(), pt3->GetE());
	//		//fMCE3 =  pt3->GetE();
	//		//fMCPx3 = pt3->GetMomentum().Px();
	//		//fMCPy3 = pt3->GetMomentum().Py();
	//		//fMCPz3 = pt3->GetMomentum().Pz();
	//		//fMCVx3 = pt3->GetStartVertex().X();
	//		//fMCVy3 = pt3->GetStartVertex().Y();
	//		//fMCVz3 = pt3->GetStartVertex().Z();
	//		//fMCT3 =  pt3->GetStartTime();
	//		//cout<<"p3 #"<<fMCPx3<<", "<<fMCPy3<<", "<<fMCPz3<<", "<<fMCE3<<std::endl;
	//		Mcp4_3.Print();
	//		fm12 = (p4_1 + p4_2).M();
	//		fm23 = (p4_2 + p4_3).M();
	//		fm31 = (p4_3 + p4_1).M();
	//		fm123 = (p4_1 + p4_2 + p4_3).M();
	//		fMCm12 = (Mcp4_1 + Mcp4_2).M();
	//		fMCm23 = (Mcp4_2 + Mcp4_3).M();
	//		fMCm31 = (Mcp4_3 + Mcp4_1).M();
	//		fMCm123 = (Mcp4_1 + Mcp4_2 + Mcp4_3).M();

	//		fMCJpsi = Mcp4_1 + Mcp4_2 + Mcp4_3;
	//		fJpsi = p4_1 + p4_2 + p4_3;
	//		tResult->Fill();
	//	}
		//std::cout<<"MC tracks #"<<fMCTrackArray->GetEntriesFast()<<std::endl;
		//for(Int_t i=0;i<fMCTrackArray->GetEntriesFast();++i){
		//	PndMCTrack* pt=((PndMCTrack*)fMCTrackArray->At(i));
		//if(pt->IsGeneratorCreated()){
		//}
		//		pt->Print(i);
		//std::cout<<"pdg #"<<pt->GetPdgCode()<<std::endl;
		//std::cout<<"Mother ID #"<<pt->GetMotherID()<<std::endl;
		//std::cout<<"IsGeneratorCreated #"<<pt->IsGeneratorCreated()<<std::endl;
		//std::cout<<"IsGeneratorDecayed #"<<pt->IsGeneratorDecayed()<<std::endl;
		//std::cout<<"Momentum #"<<pt->GetMomentum()<<std::endl;
		//std::cout<<"Vertex #"<<pt->GetStartVertex()<<std::endl;
		//std::cout<<"StartTime#"<<pt->GetStartTime()<<std::endl;
		//	}

	//
  DetectorId detIds[15] = { kDCH,kDRC,kDSK,kEMC,kGEM,kLUMI,kMDT,kMVD,kRPC,kSTT,kTPC,kTOF,kFTS,kHYPG,kHYP };
	PndMCTrack* pt1 =((PndMCTrack*)fMCTrackArray->At(0));
	Mcp4_1 = TLorentzVector(pt1->GetMomentum(), pt1->GetE());
	fTheta1 = Mcp4_1.Theta();//TMath::ACos(Mcp4_1[2]/pt1->GetMomentum().Mag());
	fPhi1 = Mcp4_1.Phi();//TMath::ATan(Mcp4_1[1]/Mcp4_1[0]);
	fHitEmc1 = pt1->GetNPoints(kEMC);
	fDet1 = -1;
	for(Int_t i=0;i<15;++i){
		//std::cout<<"detIds["<<i<<"]="<<detIds[i]<<std::endl;
		if(detIds[i] == kEMC){
			if(pt1->GetNPoints(detIds[i])>0)
			{fDet1 = detIds[i];break;}
		}else{
			if(pt1->GetNPoints(detIds[i])>0)
			{fDet1 = detIds[i];break;}
		}
	}

	PndMCTrack* pt2 =((PndMCTrack*)fMCTrackArray->At(1));
	Mcp4_2 = TLorentzVector(pt2->GetMomentum(), pt2->GetE());
	fTheta2 = Mcp4_2.Theta();//TMath::ACos(Mcp4_2[2]/pt2->GetMomentum().Mag());
	fPhi2 = Mcp4_2.Phi();//TMath::ATan(Mcp4_2[1]/Mcp4_2[0]);
	fHitEmc2 = pt2->GetNPoints(kEMC);
	fDet2 = -1;
	for(Int_t i=0;i<15;++i){
		//std::cout<<"detIds["<<i<<"]="<<detIds[i]<<std::endl;
		if(detIds[i] == kEMC){
			if(pt2->GetNPoints(detIds[i])>0)
			{fDet2 = detIds[i];break;}
		}else{
			if(pt2->GetNPoints(detIds[i])>0)
			{fDet2 = detIds[i];break;}
		}
	}

	PndMCTrack* pt3 =((PndMCTrack*)fMCTrackArray->At(2));
	Mcp4_3 = TLorentzVector(pt3->GetMomentum(), pt3->GetE());
	fTheta3 = Mcp4_3.Theta();//TMath::ACos(Mcp4_3[2]/pt3->GetMomentum().Mag());
	fPhi3 = Mcp4_3.Phi();//TMath::ATan(Mcp4_3[1]/Mcp4_3[0]);
	fHitEmc3 = pt3->GetNPoints(kEMC);
	fDet3 = -1;
	for(Int_t i=0;i<15;++i){
		//std::cout<<"detIds["<<i<<"]="<<detIds[i]<<std::endl;
		if(detIds[i] == kEMC){
			if(pt3->GetNPoints(detIds[i])>0)
			{fDet3 = detIds[i];break;}
		}else{
			if(pt3->GetNPoints(detIds[i])>0)
			{fDet3 = detIds[i];break;}
		}
	}

	tResult->Fill();

	//std::cout<<"trk 0#";
	const Double_t BarrelInner = 22./180.*TMath::Pi();
	const Double_t BarrelOuter = 140./180.*TMath::Pi();
	const Double_t ForwardInner = 5.2/180.*TMath::Pi();
	const Double_t ForwardOuter = 23.2/180.*TMath::Pi();
	const Double_t BackwardInner = 151.4/180.*TMath::Pi();
	const Double_t BackwardOuter = 169.7/180.*TMath::Pi();

	if( fTheta1 > 0. && fTheta1 < ForwardInner)//shashylik
		{
		std::cout<<"go to shashylik "<<std::endl;

		}else if( fTheta1 > ForwardInner && fTheta1 < ForwardOuter)
		{
		std::cout<<"go to forward "<<std::endl;
		std::cout<<"MVD#"<<	pt1->GetNPoints(kMVD)<<", STT#"<<pt1->GetNPoints(kSTT)
		<<", GEM#"<<pt1->GetNPoints(kGEM)<<", EMC#"<<pt1->GetNPoints(kEMC)<<std::endl;
		}else if( fTheta1 > BarrelInner && fTheta1 < BarrelOuter)
		{
		std::cout<<"go to barrel "<<std::endl;
		std::cout<<"MVD#"<<	pt1->GetNPoints(kMVD)<<", STT#"<<pt1->GetNPoints(kSTT)
		<<", GEM#"<<pt1->GetNPoints(kDCH)<<", EMC#"<<pt1->GetNPoints(kEMC)<<std::endl;
		}else if(fTheta1 > BackwardInner && fTheta1 < BackwardOuter)
		{
		std::cout<<"go to backward "<<std::endl;
		std::cout<<"MVD#"<<	pt1->GetNPoints(kMVD)<<", STT#"<<pt1->GetNPoints(kSTT)
		<<", EMC#"<<pt1->GetNPoints(kEMC)<<std::endl;
		}else{
		std::cout<<"go to empty"<<std::endl;
		};

	if (fVerbose>2){
		timer.Stop();
		Double_t rtime = timer.RealTime();
		Double_t ctime = timer.CpuTime();
		cout << "PndEmcMcTruthWriter, Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	}
	//fRecoHitArray->Delete();
	evtNo ++;
}

void PndEmcMcTruthWriter::SetParContainers() {

	// Get run and runtime database
	FairRun* run = FairRun::Instance();
	if ( ! run ) Fatal("SetParContainers", "No analysis run");

	FairRuntimeDb* db = run->GetRuntimeDb();
	if ( ! db ) Fatal("SetParContainers", "No runtime database");

	// Get Emc digitisation parameter container
	//fDigiPar = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar");

	// Get Emc reconstruction parameter container
	//fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar");

}

void PndEmcMcTruthWriter::SetStorageOfData(Bool_t val)
{
	fStoreHist= val;
	return;
}
void PndEmcMcTruthWriter::FinishTask()
{
	TFile* oldFile = gFile;
	gFile = fResult;
	gFile->cd();
	tResult->Write();
	fResult->Close();
	gFile = oldFile;
}

ClassImp(PndEmcMcTruthWriter)

	/*for(Int_t i=0;i<goodHits.size(); ++i){
		PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(goodHits[i].iHit);
		TVector3 hit_mom = theHit->where();

		PndMCTrack* pt1 =((PndMCTrack*)fMCTrackArray->At(0));
		PndMCTrack* pt2 =((PndMCTrack*)fMCTrackArray->At(1));
		PndMCTrack* pt3 =((PndMCTrack*)fMCTrackArray->At(2));
	//Double_t angle1 = fabs(theHit->GetEnergyCorrected() - pt1->GetE());
	//Double_t angle2 = fabs(theHit->GetEnergyCorrected() - pt2->GetE());
	//Double_t angle3 = fabs(theHit->GetEnergyCorrected() - pt3->GetE());
	Double_t angle1 = hit_mom.Angle(pt1->GetMomentum());
	Double_t angle2 = hit_mom.Angle(pt2->GetMomentum());
	Double_t angle3 = hit_mom.Angle(pt3->GetMomentum());
	//std::cout<<"angle1 #"<<angle1<<", angle2 #"<<angle2 <<", angle3 #"<<angle3<<std::endl;
	if(i == 0){
	fEnergy1 = theHit->energy();
	fEnergy1C = theHit->GetEnergyCorrected();
	v_1 = theHit->where();
	TVector3 pos = theHit->position();
	Double_t theta = v_1.Theta();
	Double_t phi = v_1.Phi();
	fTheta1 = v_1.Theta();
	fPhi1 = v_1.Phi();

	p4_1 = TLorentzVector(fEnergy1C*TMath::Sin(theta)*TMath::Cos(phi)
	, fEnergy1C*TMath::Sin(theta)*TMath::Sin(phi)
	, fEnergy1C*TMath::Cos(theta), fEnergy1C);
	//fX1 = theHit->x();
	//fY1 = theHit->y();
	//fZ1 = theHit->z();
	fZ201 = theHit->Z20();
	fZ531 = theHit->Z53();
	fLat1 = theHit->LatMom();
	fModule1 = theHit->GetModule();
	}
	if(i == 1){
	fEnergy2 = theHit->energy();
	fEnergy2C = theHit->GetEnergyCorrected();
	v_2 = theHit->where();
	Double_t theta = v_2.Theta();
	Double_t phi = v_2.Phi();
	fTheta2 = v_2.Theta();
	fPhi2 = v_2.Phi();
	p4_2 = TLorentzVector(fEnergy2C*TMath::Sin(theta)*TMath::Cos(phi)
	, fEnergy2C*TMath::Sin(theta)*TMath::Sin(phi)
	, fEnergy2C*TMath::Cos(theta), fEnergy2C);
	//fX2 = theHit->x();
	//fY2 = theHit->y();
	//fZ2 = theHit->z();
	fZ202 = theHit->Z20();
	fZ532 = theHit->Z53();
	fLat2 = theHit->LatMom();
	fModule2 = theHit->GetModule();
	}
	if(i == 2){
	fEnergy3 = theHit->energy();
	fEnergy3C = theHit->GetEnergyCorrected();
	v_3 = theHit->where();
	Double_t theta = v_3.Theta();
	Double_t phi = v_3.Phi();
	fTheta3 = v_3.Theta();
	fPhi3 = v_3.Phi();
	p4_3 = TLorentzVector(fEnergy3C*TMath::Sin(theta)*TMath::Cos(phi)
	, fEnergy3C*TMath::Sin(theta)*TMath::Sin(phi)
	, fEnergy3C*TMath::Cos(theta), fEnergy3C);
	//fX3 = theHit->x();
	//fY3 = theHit->y();
	//fZ3 = theHit->z();
	fZ203 = theHit->Z20();
	fZ533 = theHit->Z53();
	fLat3 = theHit->LatMom();
	fModule3 = theHit->GetModule();
	}

if(i >= 3){
	//PndEmcBump* theHit = (PndEmcBump*) fRecoHitArray->At(goodHits[i].iHit);
	fOtherEnergy += theHit->energy();
	fOtherEnergyC += theHit->GetEnergyCorrected();
}
}*/