// $Id: hgeantkineana.cc,v 1.13 2009-07-15 11:39:21 halo Exp $
// Last update by Thomas Eberl: 03/04/29 15:36:38
//
using namespace std;
#include "hgeantkineana.h"
#include "hgeantrich.h"
#include "htrackinfo.h"
#include "hgeantkine.h"
#include "hhitmatchsim.h"
#include "hdihitmatch.h"
#include "hconstant.h"
#include "hspectrometer.h"
#include "hdetector.h"
#include "hrichdetector.h"
#include "hparset.h"
#include "hrichhitsim.h"
#include "hrichutilfunc.h"
#include "hruntimedb.h"
#include "hevent.h"
#include "hcategory.h"
#include "hiterator.h"
#include "hdebug.h"
#include "hades.h"
#include "richdef.h"
#include "hlinearcategory.h"
#include "hmatrixcategory.h"
#include "hgeomvector2.h"
#include "hrichgeometrypar.h"
#include "hrichcut.h"
#include <fstream>
#include <stdlib.h>
//_HADES_CLASS_DESCRIPTION 
///////////////////////////////////////////////////////////
// HGeantKineAna
//
// this class analyses HGeantKine
// and fills histograms
///////////////////////////////////////////////////////////

ClassImp(HGeantKineAna)

HGeantKineAna::HGeantKineAna(const Text_t *name,const Text_t *title) :
  HReconstructor(name,title)
{
     
}
HGeantKineAna::HGeantKineAna(const Text_t *name,const Text_t *title,const  Char_t* filename,const Char_t* evtgen) :
  HReconstructor(name,title)
{
    pFileName  = filename;
    pEvtGen    = evtgen;
}
HGeantKineAna::HGeantKineAna()
{

}


HGeantKineAna::~HGeantKineAna(void) 
{

}

Bool_t HGeantKineAna::init() {
    if (gHades) {
	HEvent *event=gHades->getCurrentEvent();
	HRuntimeDb *rtdb=gHades->getRuntimeDb();
	HSpectrometer *spec=gHades->getSetup();
	HRichDetector *rich = (HRichDetector*)spec->getDetector("Rich");
	if (event && rtdb && rich) {
	    	    
// 	    HRichGeometryPar *pGeoParl =
// 		(HRichGeometryPar*)rtdb->getContainer("RichGeometryParameters");
// 	    setGeoPar(pGeoParl);
//             if (!pGeoPar) Error("init","no geometry cntr defined");
	    
	    fRichPID=gHades->getCurrentEvent()->getCategory(catRichHit);
	    if (!fRichPID) Warning("init","no HRichHit cat");
    
	    if (fRichPID) fRichIter = (HIterator*) fRichPID->MakeIterator();
	    if(!fRichIter) Warning("init","no RichHitSim iter defined");

	    // HGEANT RICH MIRROR INFO
	    fGeantRichMirrorCat = (HMatrixCategory*)(gHades->getCurrentEvent()
                             ->getCategory(catRichGeantRaw+2));  
	    if (!fGeantRichMirrorCat)
	    {
		Error("HRichCorrelatorSim::init",
		      "no GEANT RICH MIRROR category available");
	    }
	    iter_mirror = (HIterator *)fGeantRichMirrorCat->MakeIterator("native");
	    // HGEANT MDC RAW INFO
	    fGeantMdcCat = (HMatrixCategory*)(gHades->getCurrentEvent()
                             ->getCategory(catMdcGeantRaw));  
	    if (!fGeantMdcCat)
	    {
		Error("HRichCorrelatorSim::init",
		      "no GEANT MDC RAW category available");
	    }
	    iter_mdcgeant = (HIterator *)fGeantMdcCat->MakeIterator("native");

	    // HGEANT KINE INFO
	    fGeantKineCat =(HLinearCategory*) event->getCategory(catGeantKine);
	    if (!fGeantKineCat) 
	    { 
		Error("HGeantKineAna::init",
		      "no GEANT Kine category available");
	    }  
	    setGeantKineCat(fGeantKineCat);
	    //	    iter_kine = (HIterator *)fGeantKineCat ->MakeIterator("native");
	    iter_kine = (HIterator *) getGeantKineCat()->MakeIterator("native");

// 	    // TRACKS
// 	    pHitMatchCat=event->getCategory(catMatchHit);
// 	    if (!pHitMatchCat) {
// 		Warning("init","No HIT MATCH category defined");
// 		return kFALSE;
// 	    }
// 	    if (pHitMatchCat) pIterMatchHit = (HIterator*)getHitMatchCat()->MakeIterator("native");
// 	    // PAIRS
//   	    pDiHitMatchCat=event->getCategory(catDiMatchHit);
//   	    if (!pDiHitMatchCat) {
//   		Warning("init","No DI HIT MATCH category defined");
//   		return kFALSE;
//   	    }
// 	    if (pDiHitMatchCat) pIterDiHitMatch = (HIterator*)getDiHitMatchCat()->MakeIterator("native");
	    
	  
	} else {
	    Error ("init","! (event && rtdb)");
	    return kFALSE; //! (event && rtdb)
	}
	
//  	iniCuts();
  	iniSwitches();
	iniControlHistos();
  	iniCounters();
	pFileOut = new TFile(pFileName.Data(),"RECREATE"); // histo file
	el = new TEventList("rings","rings",100);
	
	return kTRUE;
    } else {
	Error ("init","! (gHades)");
	return kFALSE; //! (gHades)		
    }
}
void HGeantKineAna::iniSwitches()
{
    kDumpIt = kTRUE; //activate dumping to stdout of kine
}
void HGeantKineAna::iniCounters()
{
    //nCounterProcessedNbEvents=0;
    nEvtsProcessed=0;
    nbRings=0;
    nbLeptons=0;
    nbDalitzPairs=0;
    nbDalitzSingle=0;
    nbTracks=0;
    nbTrackPairs=0;
    nbgMatchedPairs=0;
    nbgMatchedTracks=0;
    nbgMatchedGPairs=0;
    nbgMatchedGTracks=0;
}

void HGeantKineAna::iniControlHistos()
{
    pHistArray = new TObjArray(12);
    med_conv = new TH1D("Medium_conv","Medium_conv",41,0,42);
    pHistArray->Add(med_conv);
    mech_conv = new TH1D("Mechanism_conv","Mechanism_conv",36,0,35);
    pHistArray->Add(mech_conv);
    par_conv = new TH1D("Parent_conv","Parent_conv",56,0,55);
    pHistArray->Add(par_conv);
    par_med_conv = new TH2D("parent_medium_conv","parent_medium_conv",56,0,55,46,0,45);
    pHistArray->Add(par_med_conv);
    par_mech_conv = new TH2D("parent_mechanism_conv","parent_mechanism_conv",56,0,55,36,0,35);
    pHistArray->Add(par_mech_conv);
    mech_mom_conv = new TH2D("mechanism_momentum_conv","mechanism_momentum_conv",36,0,35,100,0,1000);
    pHistArray->Add(mech_mom_conv);
    g_conv_opangle = new TH1F("pi0_2gamma_opangle","pi0->2*gamma_opangle",180,0,180);
    pHistArray->Add(g_conv_opangle);
    deltaP_opangle = new TH2F("delta_AbsP_opangle","delta_Abs(P)__opangle",100,0,1000,180,0,180);
    pHistArray->Add(deltaP_opangle);
    // *********************************************************
    med_dalitz = new TH1D("Medium_dalitz","Medium_dalitz",41,0,42);
    pHistArray->Add(med_dalitz);
    mech_dalitz = new TH1D("Mechanism_dalitz","Mechanism_dalitz",36,0,35);
    pHistArray->Add(mech_dalitz);
    par_dalitz = new TH1D("Parent_dalitz","Parent_dalitz",56,0,55);
    pHistArray->Add(par_dalitz);
    par_med_dalitz = new TH2D("parent_medium_dalitz","parent_medium_dalitz",56,0,55,46,0,45);
    pHistArray->Add(par_med_dalitz);
    par_mech_dalitz = new TH2D("parent_mechanism_dalitz","parent_mechanism_dalitz",56,0,55,36,0,35);
    pHistArray->Add(par_mech_dalitz);
    mech_mom_dalitz = new TH2D("mechanism_momentum_dalitz","mechanism_momentum_dalitz",36,0,35,100,0,1000);
    pHistArray->Add(mech_mom_dalitz);
    // ***********************************************************
    // *** Kine histograms ***
    hidkine = new TH1D("ID_kine","ID_kine",50,0,50);
    pHistArray->Add(hidkine);
    hmedkine = new TH1D("Medium_kine","Medium_kine",41,0,42);
    pHistArray->Add(hmedkine);
    hmechkine = new TH1D("Mechanism_kine","Mechanism_kine",36,0,35);
    pHistArray->Add(hmechkine);
    hmomkine = new TH1D("Momentum_kine","Momentum_kine",25,0,500);
    pHistArray->Add(hmomkine);
    hmomthetakine = new TH2D("Momentum_theta_kine","Momentum_theta_kine",180,0,180,100,0,500);
    pHistArray->Add(hmomthetakine);
    hmompairkine = new TH2D("Momentum_pair_kine","Momentum_pair_kine",100,0,500,100,0,500);
    pHistArray->Add(hmompairkine);
    hthetapairkine = new TH2D("polangle_pair_kine","polangle_pair_kine",180,0,180,180,0,180);
    pHistArray->Add(hthetapairkine);
    hvertexkine = new TH3D("Vertex_kine","Vertex_kine",100,-10,10,
			   100,-10,10,100,-10,10);
    pHistArray->Add(hvertexkine);
    hparentkine  = new TH1D("parent_kine","parent_kine",50,0,50);
    pHistArray->Add(hparentkine);
    // tmp histo to check dist of center
    hcntdistx = new TH1D("xdist","xdist",21,-10,10);
    pHistArray->Add(hcntdistx);
    hcntdisty = new TH1D("ydist","ydist",21,-10,10);
    pHistArray->Add(hcntdisty);
    hcntdistxy = new TH2D("xydist","xydist",21,-10,10,21,-10,10);
    pHistArray->Add(hcntdistxy);
    hweight = new TH1D("weight","weight",40,0,1.1);
    pHistArray->Add(hweight);
    hweightdist = new TH2D("weightdist","weightdist",40,0,1.1,22,-1,10.);
    pHistArray->Add(hweightdist);

    hpi0invmass = new TH1D("pi0mass","pi0mass",300,0,300);
    pHistArray->Add(hpi0invmass);
    hdalitzopang = new TH1D("Dalitzopang","Dalitzopang",101,0,100);
    pHistArray->Add(hdalitzopang);
    hdalitzinvmass = new TH1D("Dalitzinvmass","Dalitzinvmass",151,0,150);
    pHistArray->Add(hdalitzinvmass);

    hgammathetatheta = new TH2D("pi0gammatheta_theta","pi0gammatheta_theta",
				180,0,180,180,0,180);
    pHistArray->Add(hgammathetatheta);
    hgammamommom = new TH2D("pi0gammamom_mom","pi0gammamom_mom",
				100,0,500,100,0,500);
    pHistArray->Add(hgammamommom);
    hgammamom1theta2 = new TH2D("pi0gammamom1_theta2","pi0gammamom1_theta2",
				100,0,500,180,0,180);
    pHistArray->Add(hgammamom1theta2);
    hgammamomtheta = new TH2D("pi0gammamom_theta","pi0gammamom_theta",
				100,0,500,180,0,180);
    pHistArray->Add(hgammamomtheta);

    // histos for efficiency
    hmomthetaallsingle = new TH2D("Momentum_theta_allsingle","Momentum_theta_allsingle",30,0,300,20,0,90);
    pHistArray->Add(hmomthetaallsingle);
    hmomthetareconsingle = new TH2D("Momentum_theta_reconsingle","Momentum_theta_reconsingle",30,0,300,20,0,90);
    pHistArray->Add(hmomthetareconsingle);

    hmomthetaallpair = new TH2D("Momentum_theta_allpair","Momentum_theta_allpair",30,0,300,20,0,90);
    pHistArray->Add(hmomthetaallpair);
    hmomthetareconpair = new TH2D("Momentum_theta_reconpair","Momentum_theta_reconpair",30,0,300,20,0,90);
    pHistArray->Add(hmomthetareconpair);

    // conv lep histos
    hpolpol = new TH2D("diff_polang_gamma_lepton1.2","diff_polang_gamma_lepton1.2",200,0,20,200,0,20);
    pHistArray->Add(hpolpol);
    hpoldiffinvmass = new TH2D("diff_polangdiff_invmass","diff_polangdiff_invmass",600,0,300,500,0,50);
    pHistArray->Add(hpoldiffinvmass);
    hmomdiffpoldiff = new TH2D("diff_polangdiff_mom","diff_polangdiff_mom",1000,0,100,500,0,50);
    pHistArray->Add(hmomdiffpoldiff);
}

Bool_t HGeantKineAna::reconPionFromGammaConversionLeptons(TObjArray* t)
{
    Bool_t isReconPi0 = kFALSE;
    if (t->GetLast()+1<4) return isReconPi0;
    for (Int_t i=0;i<(t->GetLast()+1);i+=4)
    {
	//cout<<"this is i "<<i<<" <  "<< t->GetLast()+1 <<endl;
	
	// use two consec. leptons and form pair
	HGeantKine *e1 = (HGeantKine*)(*t)[i];
	//HRichUtilFunc::dumpKine(e1);
	HGeantKine *e2 = 0;
	if (i+1 <= t->GetLast()) e2 = (HGeantKine*)(*t)[i+1];
	else return isReconPi0;
	Int_t aTrack1, aID1;
	Int_t aPar1,aMed1,aMech1;
	Int_t aTrack2, aID2;
	Int_t aPar2,aMed2,aMech2;

	if (e1)
	{
	    e1->getParticle(aTrack1,aID1);
	    e1->getCreator(aPar1,aMed1,aMech1);
	}

	if (e2)
	{
	    e2->getParticle(aTrack2,aID2);
	    e2->getCreator(aPar2,aMed2,aMech2);
	}

	Float_t p1sumx,p1sumy,p1sumz;
	p1sumx=p1sumy=p1sumz=0;
	Float_t t1mean,p1mean;
	t1mean=p1mean=0;
	if ( e1 && e2 &&
	     aID1+aID2 == 5 && 
	     aPar1==aPar2 && 
	     HRichUtilFunc::getParentParID(aPar1,(HLinearCategory*)fGeantKineCat)==7 && // it was an e+/e- comb.
	     HRichUtilFunc::calcOpeningAngleKine(e1,e2)>0.
	     )
	{
	    Float_t apx1, apy1, apz1;
	    e1->getMomentum(apx1,apy1,apz1);
	    Float_t apx2, apy2, apz2;
	    e2->getMomentum(apx2,apy2,apz2);
	    p1sumx=apx1+apx2;
	    p1sumy=apy1+apy2;
	    p1sumz=apz1+apz2;
 	    Float_t theta1,phi1;
 	    HRichUtilFunc::calcParticleAnglesKine(e1,theta1,phi1);
 	    Float_t theta2,phi2;
 	    HRichUtilFunc::calcParticleAnglesKine(e2,theta2,phi2);
	    t1mean = (theta1+theta2)/2.;
	    p1mean = (phi1+phi2)/2.;
	    HGeantKine *gamma = HRichUtilFunc::findParent(e2,(HLinearCategory*)fGeantKineCat);
	    Float_t gammat,gammap;
 	    HRichUtilFunc::calcParticleAnglesKine(gamma,gammat,gammap);
	    gamma->getMomentum(apx1,apy1,apz1);
	    hpolpol->Fill(TMath::Abs(gammat-theta1),(TMath::Abs(gammat-theta2)));
// 	    cout << "theta diff arithm. mean e+e- and gamma: " << gammat-t1mean <<endl;
// 	    cout << "phi   diff arithm. mean e+e- and gamma: " << gammap-p1mean << endl;
// 	    cout << apx1-p1sumx << "  " << apy1-p1sumy << "  " << apz1-p1sumz << endl;

	}

	HGeantKine *e3 = 0;
	HGeantKine *e4 = 0;
	if (i+2 <= t->GetLast()) e3 = (HGeantKine*)(*t)[i+2];
	else return isReconPi0;
	if (i+3 <= t->GetLast()) e4 = (HGeantKine*)(*t)[i+3];
	else return isReconPi0;
	Int_t aTrack3, aID3;
	Int_t aPar3,aMed3,aMech3;
	Int_t aTrack4, aID4;
	Int_t aPar4,aMed4,aMech4;
	if (e3)
	{
	    e3->getParticle(aTrack3,aID3);
	    e3->getCreator(aPar3,aMed3,aMech3);
	}
	if (e4)
	{
	    e4->getParticle(aTrack4,aID4);
	    e4->getCreator(aPar4,aMed4,aMech4);
	}

	Float_t p2sumx,p2sumy,p2sumz;
	p2sumx=p2sumy=p2sumz=0;
	Float_t t2mean,p2mean;
	t2mean=p2mean=0;
	if ( e3 && e4 && 
	     aID3+aID4 == 5 && 
	     aPar3==aPar4 && 
	     HRichUtilFunc::getParentParID(aPar3,(HLinearCategory*)fGeantKineCat)==7 && // it was an e+/e- comb.
	     HRichUtilFunc::calcOpeningAngleKine(e3,e4)>0.
	     )
	{
	    Float_t apx1, apy1, apz1;
	    e3->getMomentum(apx1,apy1,apz1);
	    Float_t apx2, apy2, apz2;
	    e4->getMomentum(apx2,apy2,apz2);
	    p2sumx=apx1+apx2;
	    p2sumy=apy1+apy2;
	    p2sumz=apz1+apz2;
 	    Float_t theta1,phi1;
 	    HRichUtilFunc::calcParticleAnglesKine(e3,theta1,phi1);
 	    Float_t theta2,phi2;
 	    HRichUtilFunc::calcParticleAnglesKine(e4,theta2,phi2);
	    t2mean = (theta1+theta2)/2.;
	    p2mean = (phi1+phi2)/2.;
	    HGeantKine *gamma = HRichUtilFunc::findParent(e3,(HLinearCategory*)fGeantKineCat);
	    Float_t gammat,gammap;
 	    HRichUtilFunc::calcParticleAnglesKine(gamma,gammat,gammap);
	    hpolpol->Fill(TMath::Abs(gammat-theta1),(TMath::Abs(gammat-theta2)));
	  //   cout << "theta diff arithm. mean e+e- and gamma: " << gammat-t2mean << "   "
// 		 << "phi   diff arithm. mean e+e- and gamma: " << gammap-p2mean << endl;
	}
	if (e1 && e2 && e3 && e4)
	{
	    Float_t p1sumtot = TMath::Sqrt(p1sumx*p1sumx+p1sumy*p1sumy+p1sumz*p1sumz); 
	    Float_t p2sumtot = TMath::Sqrt(p2sumx*p2sumx+p2sumy*p2sumy+p2sumz*p2sumz); 
	    Float_t opang    = HRichUtilFunc::calcOpeningAngleT(t1mean,p1mean,t2mean,p2mean);
	    Float_t invmass  = 2.*sin(opang/57.3/2.)*sqrt(p1sumtot*p2sumtot);
	    //	cout<<"Invariant mass calculated: "<<invmass<<endl;
	    if (invmass>0.) 
	    {

		//hpi0invmass->Fill(invmass);
		Float_t theta1,phi1;
		HRichUtilFunc::calcParticleAnglesKine(e1,theta1,phi1);
		Float_t theta2,phi2;
		HRichUtilFunc::calcParticleAnglesKine(e2,theta2,phi2);
		Float_t theta3,phi3;
		HRichUtilFunc::calcParticleAnglesKine(e3,theta3,phi3);
		Float_t theta4,phi4;
		HRichUtilFunc::calcParticleAnglesKine(e4,theta4,phi4);
		Float_t mom1,mom2,mom3,mom4;
		mom1 = e1->getTotalMomentum();
		mom2 = e2->getTotalMomentum();
		mom3 = e3->getTotalMomentum();
		mom4 = e4->getTotalMomentum();
		HGeantKine *gamma1 = HRichUtilFunc::findParent(e1,(HLinearCategory*)fGeantKineCat);
		HGeantKine *gamma2 = HRichUtilFunc::findParent(e3,(HLinearCategory*)fGeantKineCat);
		Float_t gamma1t,gamma1p,gamma2t,gamma2p;
		HRichUtilFunc::calcParticleAnglesKine(gamma1,gamma1t,gamma1p);
		HRichUtilFunc::calcParticleAnglesKine(gamma2,gamma2t,gamma2p);
		hpoldiffinvmass->Fill(invmass,fabs(gamma1t-theta1)-fabs(gamma1t-theta2));
		hpoldiffinvmass->Fill(invmass,fabs(gamma2t-theta3)-fabs(gamma2t-theta4));

		hmomdiffpoldiff->Fill(fabs(mom1-mom2),fabs(gamma1t-theta1)-fabs(gamma1t-theta2));
		hmomdiffpoldiff->Fill(fabs(mom3-mom4),fabs(gamma2t-theta3)-fabs(gamma2t-theta4));
	    }
	    if (invmass>130. && invmass<140.) 
	    {
		// HRichUtilFunc::dumpKine(e1);
// 		HRichUtilFunc::dumpKine(e2);
// 		HRichUtilFunc::dumpKine(e3);
// 		HRichUtilFunc::dumpKine(e4);
// 		cout<<"^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^"<<endl;
		isReconPi0=kTRUE;

	    }
	}
    }


    return isReconPi0;
}

TObjArray* HGeantKineAna::getPionGammaConversionLeptons()
{

    TObjArray *leps = new TObjArray(4);
    // find pi0 

    HGeantKine * pion =0;
    iter_kine->Reset();

    while((pion=(HGeantKine *)iter_kine->Next())!=0)
    {
	Int_t aID, aTrack;
	pion->getParticle(aTrack,aID);
	if (aID==7) // is pi0
	{
	    //cout<<"is pion"<<endl;
	    // find gammas 
	    HIterator* iter_gamma = (HIterator*)(fGeantKineCat->MakeIterator("native"));
	    iter_gamma->Reset();
	    HGeantKine *gamma=0;
	    while((gamma=(HGeantKine *)iter_gamma->Next())!=0)
	    {
		
		Int_t aTrackGamma, aIDGamma;
		Int_t aParGamma,aMedGamma,aMechGamma;
		gamma->getParticle(aTrackGamma,aIDGamma);
		gamma->getCreator(aParGamma,aMedGamma,aMechGamma);
		if (aIDGamma==1 && aMechGamma==5 && aParGamma==aTrack)
		{
		    //cout<<"is gamma from pion"<<endl;
		    // find e+e- for gamma
		    HGeantKine *gamma1 = HRichUtilFunc::getSecondPionDecayGamma(gamma,(HLinearCategory*)fGeantKineCat);
		    if (gamma1)
		    {
			Float_t mom1,mom2;
			mom1 = gamma->getTotalMomentum();
			mom2 = gamma1->getTotalMomentum();
			Float_t opang = HRichUtilFunc::calcOpeningAngleKine(gamma,gamma1);
			Float_t invmass  = 2.*sin(opang/57.3/2.)*sqrt(mom1*mom2);
			hpi0invmass->Fill(invmass);
		    }
		    HIterator* iter_ele = (HIterator*)(fGeantKineCat->MakeIterator("native"));
		    iter_ele->Reset();
		    HGeantKine *ele=0;
		    while((ele=(HGeantKine *)iter_ele->Next())!=0)
		    {
			
			Int_t aTrackEle, aIDEle;
			Int_t aParEle,aMedEle,aMechEle;
			ele->getParticle(aTrackEle,aIDEle);
			ele->getCreator(aParEle,aMedEle,aMechEle);
			if ((aIDEle==2 || aIDEle==3) && aMechEle==6 && aParEle==aTrackGamma)
			{
			    // cut on acceptance:
			    //HGeantRichMirror *mirr=0;
			    //iter_mirror->Reset();
			    //while((mirr=(HGeantRichMirror *)iter_mirror->Next())!=0)
			    //{
			    //Int_t mirrTrack,aID;
			    //mirr->getTrack(mirrTrack,aID);
			    //if (mirrTrack==aTrackEle) //lepton passed mirror
			    //{
			    //    Float_t theta,phi;
			    //    HRichUtilFunc::calcParticleAnglesKine(ele,theta,phi);
				    //cout<<"is ele from gamma"<<endl;
			    //    if (
			    //	theta>18.&& theta<85. && 
			    //	ele->getTotalMomentum()>50.
			    //	) 
			    //    {
					leps->Add(ele);
					//HRichUtilFunc::dumpKine(ele);
					//    }
					//    break;
				}
			    }
			}
	    }
		    
	}
		
    }

	
	

    return leps;
}

void HGeantKineAna::findPionDalitzLeptons()
{
    HGeantKine * kine =0;
    // loop over kine container
    iter_kine->Reset();
    while((kine=(HGeantKine *)iter_kine->Next())!=0)
    {
	Int_t aTrack, aID;
	Int_t aPar, aMed, aMech;
	//Float_t ax, ay, az;
	//Float_t apx, apy, apz;
	//Float_t aInfo, aWeight;
	Float_t ptot;
	kine->getParticle(aTrack,aID);
	kine->getCreator(aPar,aMed,aMech);
	//kine->getVertex(ax,ay,az);
	//kine->getMomentum(apx,apy,apz);
	//kine->getGenerator(aInfo,aWeight);
	ptot=kine->getTotalMomentum();

	// aID==2 is e+
	// aID==3 is e-
	// aMech==6 is photo pair production
	// aMed==8 is RICH radiator ?? aMed==14 is target ??
	// aMed==9 is RICH window ?? aMed==10 is RICH carb shell ??
	// aMed==13 is RICH glass mirror ??
//  	if ( aID == 2 && aMech == 6 && 
//  	     //  (aMed == 8 || aMed == 14 || aMed == 9 
//  	     // || aMed == 10 || aMed == 13)) pPositrons->Add(kine);
//  	     aMed == 9) pPositrons->Add(kine);
//  	if ( aID == 3 && aMech == 6 && 
//  	     //(aMed == 8 || aMed == 14 || aMed == 9
//  	     // || aMed == 10 || aMed == 13)) pElectrons->Add(kine);
//  	     aMed == 9) pElectrons->Add(kine);
	//	if ( ( aID == 2 || aID==3 ) && aMech != 10 )  dumpKine(kine);


	if ( aMech==5 && (aID == 2 || aID==3) && (aMed>=8 && aMed<=19)) // direct particle decay
	{
	    HGeantKine *kine_parent1 = HRichUtilFunc::findParent(kine,(HLinearCategory*) getGeantKineCat());
	    if(kine_parent1){//parent should be pion
		Float_t theta1,phi1;
		theta1=phi1=0.;
		HRichUtilFunc::calcParticleAnglesKine(kine_parent1,theta1,phi1);
		Int_t aTrackp1, aIDp1;
		kine_parent1->getParticle(aTrackp1,aIDp1);
		//if(aIDp1==7 && (theta1>=10. && theta1<=90.))//pion
		if(aIDp1==7)//pion, but no emission angle cut
		{ //in acceptance
		    HGeantKine *secondlepton = 
			HRichUtilFunc::getSecondPionDalitzLepton(kine,(HLinearCategory*) getGeantKineCat());
		    if (secondlepton)
		    {
			Float_t deltaptot = 
			    TMath::Abs(ptot-secondlepton->getTotalMomentum());
			Float_t opangle = 
			    HRichUtilFunc::calcOpeningAngleKine(kine,secondlepton);
			deltaP_opangle->Fill(deltaptot,opangle);
		    }
		    med_dalitz->Fill(aMed);
		    mech_dalitz->Fill(aMech);
		    par_dalitz->Fill(aIDp1);
		    par_med_dalitz->Fill(aID,aMed);
		    par_mech_dalitz->Fill(aID,aMech);
		    mech_mom_dalitz->Fill(aMech,ptot);
		}//pi0
	    }//kine_parent1
	}//Mec==5 and lepton
	//	    if (kDumpIt) dumpKine(kine);
	//  if (kDumpIt) dumpKine(findParent(kine));

    } // end while kine loop


}
void HGeantKineAna::findPionGammaConversionLeptons()
{

    HGeantKine * kine =0;
    // loop over kine container
    iter_kine->Reset();
    while((kine=(HGeantKine *)iter_kine->Next())!=0)
    {
	Int_t aTrack, aID;
	Int_t aPar, aMed, aMech;
	//Float_t ax, ay, az;
	//Float_t apx, apy, apz;
	//Float_t aInfo, aWeight;
	Float_t ptot;
	kine->getParticle(aTrack,aID);
	kine->getCreator(aPar,aMed,aMech);
	//kine->getVertex(ax,ay,az);
	//kine->getMomentum(apx,apy,apz);
	//kine->getGenerator(aInfo,aWeight);
	ptot=kine->getTotalMomentum();

	// aID==2 is e+
	// aID==3 is e-
	// aMech==6 is photo pair production
	// aMed==8 is RICH radiator ?? aMed==14 is target ??
	// aMed==9 is RICH window ?? aMed==10 is RICH carb shell ??
	// aMed==13 is RICH glass mirror ??
//  	if ( aID == 2 && aMech == 6 && 
//  	     //  (aMed == 8 || aMed == 14 || aMed == 9 
//  	     // || aMed == 10 || aMed == 13)) pPositrons->Add(kine);
//  	     aMed == 9) pPositrons->Add(kine);
//  	if ( aID == 3 && aMech == 6 && 
//  	     //(aMed == 8 || aMed == 14 || aMed == 9
//  	     // || aMed == 10 || aMed == 13)) pElectrons->Add(kine);
//  	     aMed == 9) pElectrons->Add(kine);
	//	if ( ( aID == 2 || aID==3 ) && aMech != 10 )  dumpKine(kine);


	if ( aMech==6 && (aID == 2 || aID==3) && (aMed>=8 && aMed<=19)) // photon pair production in target or Rich
	{
	    HGeantKine *kine_parent1 = HRichUtilFunc::findParent(kine,(HLinearCategory*) getGeantKineCat());
	    if(kine_parent1){
		
		Int_t aTrackp1, aIDp1;
		Int_t aParp1,aMedp1,aMechp1;
		kine_parent1->getParticle(aTrackp1,aIDp1);
		kine_parent1->getCreator(aParp1,aMedp1,aMechp1);
		if(aIDp1==1 && aMechp1==5)//gamma
		{
		    
		    HGeantKine *kine_parent2 = HRichUtilFunc::findParent(kine_parent1,(HLinearCategory*) getGeantKineCat());
		    if(kine_parent2){
			Int_t aTrackp2, aIDp2;
			//  	    cout<<"before dumpKine "<<kine_parent2<<endl;
			//			  	    dumpKine(kine_parent2);
			kine_parent2->getParticle(aTrackp2,aIDp2);
			Float_t theta2,phi2;
			theta2=phi2=0.;
			HRichUtilFunc::calcParticleAnglesKine(kine_parent2,theta2,phi2);
			if(aIDp2==7 && (theta2>=10. && theta2<=90.))
			    //if(aIDp2==7)
			{//pi0 in acceptance
			    HGeantKine *secondgamma = HRichUtilFunc::getSecondPionDecayGamma(kine_parent1,(HLinearCategory*) getGeantKineCat());
			    
			    if (secondgamma) g_conv_opangle->Fill(HRichUtilFunc::calcOpeningAngleKine(kine_parent1,secondgamma));

			    med_conv->Fill(aMed);
			    mech_conv->Fill(aMech);
			    par_conv->Fill(aIDp2);
			    par_med_conv->Fill(aID,aMed);
			    par_mech_conv->Fill(aID,aMech);
			    mech_mom_conv->Fill(aMech,ptot);
			}//pi0
		    }//kine_parent2
		}//gamma
	    }//kine_parent1
	}//Mec==6 and lepton
	//	    if (kDumpIt) dumpKine(kine);
	//  if (kDumpIt) dumpKine(findParent(kine));

    } // end while kine loop


}

Int_t HGeantKineAna::execute()
{
//     cout<<" **** Event number "<<nEvtsProcessed<<" -----------------"<<endl;
    nEvtsProcessed++;
    //cout<<"in execute"<<endl;
    if ( nEvtsProcessed%5000==0) 
    {

	cout<<nEvtsProcessed
	    <<" EVTS PROCESSED ***** saving histos"<<endl;
	HRichUtilFunc::saveHistos(pFileOut,pHistArray);
    }

    TObjArray *leps = getPionGammaConversionLeptons();
    //    if (leps->GetLast()+1 > 4){
    // 	cout<<"size of array: "<<leps->GetLast()+1<<endl;

// 	for (Int_t i=0;i<(leps->GetLast()+1);i++)
// 	{
// 	    //	    HGeantKine *ele = (HGeantKine*)(*leps)[i];
// // 	    if (ele) HRichUtilFunc::dumpKine(ele);	
// 	}
	
//    }
Bool_t kWriteEvt = reconPionFromGammaConversionLeptons(leps);
    delete leps;

//      TObjArray *pElectrons = new TObjArray(25);
//      TObjArray *pPositrons = new TObjArray(25);

//      findPionGammaConversionLeptons();
//      findPionDalitzLeptons();

//     TObjArray dalitzleptons(1);
//     if(findAcceptedDalitzLeptons(dalitzleptons)) 
// 	checkAcceptedDalitzLeptons(dalitzleptons);


    //showLeptonsFromPi0Decays("dalitz");
    //    showGammasFromPi0Decays();
    // ===========================================================


    // loop over selected e- and e+ (kine objs) and fill histos 

//      Int_t aParE,aParP,aMedE,aMedP,aMechE,aMechP,aMedG,aMechG,aParG;
//      aParE=aParP=aMedE=aMedP=aMechE=aMechP=aMedG=aMechG=aParG=-1;
//      Float_t theta,phi;
//      theta=phi =-1.;
//      for (Int_t i=0;i<(pPositrons->GetLast()+1);i++)
//      {
//  	HGeantKine *positron = (HGeantKine*)(*pPositrons)[i];
//  	if (positron) 
//  	{
//  	    positron->getCreator(aParP,aMedP,aMechP);
//  	    for (Int_t j=0;j<(pElectrons->GetLast()+1);j++)
//  	    {
//  		HGeantKine *electron = (HGeantKine*)(*pElectrons)[j];
//  		if (electron)  
//  		{
//  		    electron->getCreator(aParE,aMedE,aMechE);
		
//  		    if (aParE==aParP && aMedE==aMedP && aMechE==aMechP)
//  		    { 
//  			//cout<<"!!! Conversion pair found !!!"<<endl;
//  			//dumpKine(electron);
//  			//dumpKine(positron);
//  			pHistConvPairSource->Fill(aMedP);
//  			pHistElePosMomConv->Fill(electron->getTotalMomentum(),
//  						 positron->getTotalMomentum());
//  			pHistOpenAngle->Fill(calcOpeningAngleV(electron,
//  							       positron));
//  			calcParticleAnglesV(electron,theta,phi);
//  			pHistElePosMomThetaConv->
//  			    Fill(electron->getTotalMomentum(),
//  				 positron->getTotalMomentum(),
//  				 theta);
			
		       
			
//  		    } // end pair check
//  		} // if electron
//  	    }//end for electrons
//  	}//if positron
//      } //end for positrons

//      aParE=aParP=aMedE=aMedP=aMechE=aMechP=aMedG=aMechG=aParG=-1;
//      Int_t aTrack,aID;
//      aTrack=aID=-1;
//      Float_t vx,vy,vz;
//      vx=vy=vz=-10000.;
//      theta=phi=-1.;
//      for (Int_t k=0;k<(pElectrons->GetLast()+1);k++)
//      {
//  	HGeantKine *electron = (HGeantKine*)(*pElectrons)[k];
//  	if (electron)  
//  	{
//  	    electron->getCreator(aParE,aMedE,aMechE);
//  	    electron->getVertex(vx,vy,vz);
//  	    calcParticleAnglesV(electron,theta,phi);
//  	    pHistConvEleMom->Fill(theta,electron->getTotalMomentum());
//  	    pHistConvPairVertex->Fill(vx,vy,vz);
//  	    HGeantKine *gamma = findParent(electron);
//  	    if (gamma) 
//  	    {
//  		gamma->getParticle(aTrack,aID);
//  		if (aID == 1 && aMechE == 6 && aParE == aTrack) // gamma
//  		{
//  		   calcParticleAnglesV(gamma,theta,phi);
//  		   pHistConvGammaMom->Fill(theta,gamma->getTotalMomentum()); 
//  		}
//  	    }
//  	}
//      }
//      vx=vy=vz=-10000.;
//      for (Int_t l=0;l<(pPositrons->GetLast()+1);l++)
//      {
//  	HGeantKine *positron = (HGeantKine*)(*pPositrons)[l];
//  	if (positron)  
//  	{
//  	    positron->getCreator(aParP,aMedP,aMechP);
//  	    positron->getVertex(vx,vy,vz);
//  	    calcParticleAnglesV(positron,theta,phi);
//  	    pHistConvPosMom->Fill(theta,positron->getTotalMomentum());
//  	    pHistConvPairVertex->Fill(vx,vy,vz);
//  	    HGeantKine *gamma = findParent(positron);
//  	    if (gamma) 
//  	    {
//  		gamma->getParticle(aTrack,aID);
//  		if (aID == 1 && aMechP == 6 && aParP == aTrack) // gamma
//  		{
//  		   calcParticleAnglesV(gamma,theta,phi);
//  		   pHistConvGammaMom->Fill(theta,gamma->getTotalMomentum()); 
//  		}
//  	    }
//  	}
//      }
//      nCounterProcessedNbEvents++;
    
//      delete pElectrons;
//      delete pPositrons;
    if (kWriteEvt) return 1;
    else return kSkipEvent; //skips writing out event to file


}

HRichHitSim* HGeantKineAna::checkRings(HGeantKine* kine)
{
    // check whether the kine lepton was seen as a ring
    // or at least as part of it.
    Int_t ret_val=0;
    fRichIter->Reset();
    HRichHitSim *r = 0;
    Int_t aTrk, aID;
    kine->getParticle(aTrk,aID);

    Int_t richInd[20];
    Int_t xdistArr[20];
    Int_t ydistArr[20];
    Float_t distArr[20];
    for (Int_t i=0;i<20;i++) richInd[i]=xdistArr[i]=ydistArr[i]=-1;
    for (Int_t j=0;j<20;j++) distArr[j]=-1.;
    Int_t ringCnt=0;
    while((r = (HRichHitSim *)fRichIter->Next()))
    {
	Int_t t1,t2,t3;
	t1=t2=t3=-1;
	t1=r->track1;t2=r->track2;t3=r->track3;
	if ((t1==aTrk) || (t2==aTrk) || (t3==aTrk) ) // the ring contains same track number as the lepton
	{
	    // check whether the lepton went through the mirror where the ring sits.
	        HGeantRichMirror *mirr=0;
		iter_mirror->Reset();
		HRichGeometryPar* pGeo = (HRichGeometryPar*) pGeoPar;
		Float_t fYShift = pGeo->getSectorShift(); // to account for shiftewd volumes in HGEANT ?!
		fYShift = fYShift/TMath::Cos(20./57.3);
		HRichPadTab* pPadsTab = ((HRichGeometryPar*)pGeo)->getPadsPar();
		HRichPad *pPad = 0;
		Int_t xRing=r->getRingCenterX();
		Int_t yRing=r->getRingCenterY();
		//	Int_t secRing=r->getSector();
		// iterate over hits on mirror and compare hit ring x,y with
		// center of reflected photons, calculate pad hit by photon ctr
		while((mirr=(HGeantRichMirror *)iter_mirror->Next())!=0)
		{
		    Int_t mirrTrack,aID;
		    mirr->getTrack(mirrTrack,aID);
		    if (mirrTrack==aTrk) //lepton passed mirror
		    {
			ret_val=1;
			//return 1;
			Float_t x=mirr->getYRing()/10.;//in mm of the padplane
			Float_t y=mirr->getXRing()/10. + fYShift;
			pPad=pPadsTab->getPad(x, // they are swapped, sic !
					      y );// divide by 10 to account for mm->cm
			if (!pPad) Error("checkRings","no pad found for coord");
			//			Int_t sec=mirr->getSector();
			Int_t xMirrPhotCM=pPad->getPadX();
			Int_t yMirrPhotCM=pPad->getPadY();
			Float_t distance = TMath::Sqrt(1.*(xMirrPhotCM-xRing)*(xMirrPhotCM-xRing)+
						       1.*(yMirrPhotCM-yRing)*(yMirrPhotCM-yRing));
			Int_t richind = getRichHitCat()->getIndex(r);
			richInd[ringCnt]=richind;
			distArr[ringCnt]=distance;
			xdistArr[ringCnt]=xMirrPhotCM-xRing;
			ydistArr[ringCnt]=yMirrPhotCM-yRing;
			ringCnt++;

			Float_t padsum = r->weigTrack1+r->weigTrack2+r->weigTrack3;
			Float_t weight =-1.;
			if (aTrk==r->track1) weight = ((Float_t)r->weigTrack1)/padsum;
			else if (aTrk==r->track2) weight = ((Float_t)r->weigTrack2)/padsum;
			else if (aTrk==r->track3) weight = ((Float_t)r->weigTrack3)/padsum;
			//if (distance>5. && weight==1.)
			//{
			    //el->Enter(nEvtsProcessed-1);//account for start event
			    //el->Print("all");
			    //r->dumpToStdout();
			    //}
			//cout<<weight<<endl;
//  			hweight->Fill(weight);
//  			hweightdist->Fill(weight,distance);
//    			hcntdistx->Fill(xMirrPhotCM-xRing);
//      			hcntdisty->Fill(yMirrPhotCM-yRing);
//      			hcntdistxy->Fill(xMirrPhotCM-xRing,yMirrPhotCM-yRing);

			// matching condition via distance
			    //if (TMath::Abs(xMirrPhotCM-xRing)<=2 && 
			    //TMath::Abs(yMirrPhotCM-yRing)<=2 && 
			    //sec==secRing)
			    //{
//  			    cout<<"*********************************"<<endl;
//  			    cout<<"xmirr: "<<xMirrPhotCM<<" ymirr: ";
//  			    cout<<yMirrPhotCM<<" sec: "<<sec<<endl;
//  			    cout<<"xring: "<<xRing<<" yring: "<<yRing<<" sec: "<<sec<<endl;
//  			    cout<<"*********************************"<<endl;
			    // ring position and cm of photons from this lepton agree
			    // this lepton was seen as a ring
			    //return a positive decision for this lepton
			    //return 1;

			    //}


		    }
		}
	}
    }// end while over rich hits

    // now find ring with minimum distance to photon center

    HRichHitSim *rs = 0;
    Int_t sorted[20]; for (Int_t k=0;k<20;k++) sorted[k]=-1;
    TMath::Sort(20, distArr, sorted,kFALSE); //increasing order

//      for (Int_t k=0;k<20;k++) cout<<distArr[sorted[k]]<<" ";
//      cout<<endl;

    Int_t nn=0;
    while (nn<19&&distArr[sorted[nn]]==-1) nn++;
    //    cout<<nn<<endl;
    if (richInd[sorted[nn]]!=-1)
    {
//  	cout<<nn<<" "<<distArr[sorted[nn]]<<endl;
//  	cout<<"rich index: "<<richInd[sorted[nn]]<<endl;

  	rs =(HRichHitSim*)getRichHitCat()
  	    ->getObject(richInd[sorted[nn]]);
	
	Float_t padsum = rs->weigTrack1+rs->weigTrack2+rs->weigTrack3;
	Float_t weight =-1.;
	if (aTrk==rs->track1) weight = ((Float_t)rs->weigTrack1)/padsum;
	else if (aTrk==rs->track2) weight = ((Float_t)rs->weigTrack2)/padsum;
	else if (aTrk==rs->track3) weight = ((Float_t)rs->weigTrack3)/padsum;
	//cout<<weight<<endl;
  	hweight->Fill(weight);
  	hweightdist->Fill(weight,distArr[sorted[nn]]);
  	hcntdistx->Fill(xdistArr[sorted[nn]]);
    	hcntdisty->Fill(ydistArr[sorted[nn]]);
  	hcntdistxy->Fill(xdistArr[sorted[nn]],ydistArr[sorted[nn]]);
    }
    return rs;

}
void HGeantKineAna::showGammasFromPi0Decays()
{

    HGeantKine * kine =0;
    // loop over kine container
    iter_kine->Reset();
    Int_t size = 20;
    Int_t *gammapairs = new Int_t[size];
    for (Int_t k=0;k<size;k++) gammapairs[k]=-2;
    while((kine=(HGeantKine *)iter_kine->Next())!=0)
    {
	Int_t i = getGeantKineCat()->getIndex(kine);
	HGeantKine *secondgamma=0;
  	if(HRichUtilFunc::isGammaFromPi0(kine,
  				  (HLinearCategory*)getGeantKineCat())
	)
	    //if(HRichUtilFunc::isGamma(kine))

	{
	    secondgamma = HRichUtilFunc::getSecondPionDecayGamma(kine,
				  (HLinearCategory*)getGeantKineCat());
	    //secondgamma = kine;
	    Int_t j = getGeantKineCat()->getIndex(secondgamma);
	    
	    if (HRichCut::isNew2Tuple(i,j,gammapairs,size) &&
		secondgamma
		) histoGammas(kine,secondgamma);


	}
    }
    delete [] gammapairs;
}
void HGeantKineAna::histoGammas(HGeantKine* gamma1,HGeantKine* gamma2)
{
    
    Float_t momgamma1 = gamma1->getTotalMomentum();
    Float_t momgamma2 = gamma2->getTotalMomentum();
    
    Float_t theta_gamma1,phi_gamma1;
    Float_t theta_gamma2,phi_gamma2;
    
    HRichUtilFunc::calcParticleAnglesKine(gamma1,theta_gamma1,phi_gamma1);
    HRichUtilFunc::calcParticleAnglesKine(gamma2,theta_gamma2,phi_gamma2);
    hgammathetatheta->Fill(theta_gamma1,theta_gamma2);//account for Double_t counting
    hgammamommom->Fill(momgamma1,momgamma2);//account for Double_t counting
    hgammamomtheta->Fill(momgamma1,theta_gamma1);    
    hgammamomtheta->Fill(momgamma2,theta_gamma2);    
    hgammamom1theta2->Fill(momgamma1,theta_gamma2);

}

void HGeantKineAna::showLeptonsFromPi0Decays(const Char_t *swt)
{

    TObjArray leps(4);


    // loop over all particles provided by GEANT
    // and sort out different kinds
    HGeantKine * kine =0;
    iter_kine->Reset();
    while((kine=(HGeantKine *)iter_kine->Next())!=0)
    {


	TString decaymode(swt);
  	if (!decaymode.CompareTo("dalitz") &&
	    HRichUtilFunc::isPi0DalitzLep(kine,
					  (HLinearCategory*)getGeantKineCat(),pEvtGen.Data()) )
  	{
  	    histoKine(kine);
  	    leps.Add(kine);
  	}

	if (!decaymode.CompareTo("conv") &&
	    HRichUtilFunc::isPi0ConvLep(kine,
					(HLinearCategory*)getGeantKineCat()) )
	{
	    histoKine(kine);
	    leps.Add(kine);
	}

    	if (!decaymode.CompareTo("lepinMDC") &&
	    HRichUtilFunc::isPi0DalitzLep(kine,(HLinearCategory*)getGeantKineCat(),
					  pEvtGen.Data()) &&
      	    HRichUtilFunc::isLepOnMDC(kine,(HMatrixCategory*)getGeantMDCCat())     )
    	{
    	    histoKine(kine);
    	}

  	if (!decaymode.CompareTo("leponmirror") &&
	    HRichUtilFunc::isPi0DalitzLep(kine,(HLinearCategory*)getGeantKineCat(),
					  pEvtGen.Data()) &&
    	    HRichUtilFunc::isLepOnMirror(kine,(HMatrixCategory*)getGeantMirrCat())     )
  	{
  	    histoKine(kine);
  	}


	if (!decaymode.CompareTo("dalitzlepinacc") &&
	    HRichUtilFunc::isPi0DalitzLep(kine,(HLinearCategory*)getGeantKineCat(),pEvtGen.Data()) &&
  	    HRichUtilFunc::isLepOnMirror(kine,(HMatrixCategory*)getGeantMirrCat())  &&
  	    HRichUtilFunc::isLepOnMDC(kine,(HMatrixCategory*)getGeantMDCCat())     )
  	{
	    histoKine(kine);
	    leps.Add(kine); // store lepton for further analysis
	    cout<<" **** lepton *************************************************"<<endl;
	    HRichUtilFunc::dumpKine(kine);
	    HRichHitSim* rs=0;
  	    if ((rs=checkRings(kine))) // ring with min dist to lepton is returned 
	    {
		
		cout<<" **** corresponding min dist ring ****************************"<<endl; 
		//rs->dumpToStdout();
		
		TObjArray tracks(5);
		// this function finds tracks for this ring and applies a quality cut 
		// to the track, they are provided in array tracks 
		Int_t nbTracks = checkRingCorrelation(rs,tracks);
		if (nbTracks)
		{
		    cout<<" ***** constructed tracks for this ring *********************"<<endl;
		    for (Int_t i=0;i<(tracks.GetLast()+1);i++)
		    {
			((HHitMatchSim*)tracks[i])->dumpToStdoutSim();
		    }

		}

		// if there were any tracks, try to find pairs
		// and provide them in the array pairs
		TObjArray pairs(3);
		if (nbTracks)
		{
                    Int_t nbpairs =0;
		    nbpairs = checkTrackPairs(tracks,pairs);
		}
		
		cout<<" ****  -------------------- ***************************************"<<endl;
	    }
	    else
	    {
		// this lepton was not recognized in the RICH
		// as no ring was found with photons originating from
		// this lepton !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

		//  		HRichUtilFunc::dumpKine(kine);
		//  		el->Enter(nEvtsProcessed-1);
		
	    } // end if ring check

  	} //end if decay mode
	
    }
    
    checkForPair(leps,swt);
    
    
}



Int_t HGeantKineAna::findAcceptedDalitzLeptons(TObjArray& dalitzleps)
{
    // this function fills Dalitz leptons from HGEantKine in dalitzleps 
    Int_t n=0; 
    HGeantKine * kine =0;
    iter_kine->Reset();
    while((kine=(HGeantKine *)iter_kine->Next())!=0)
    {
	//HRichUtilFunc::dumpKine(kine);
	if (HRichUtilFunc::isPi0DalitzLep(kine,(HLinearCategory*)getGeantKineCat(),pEvtGen.Data()) &&
	    HRichUtilFunc::isLepOnMirror(kine,(HMatrixCategory*)getGeantMirrCat())  &&
	    HRichUtilFunc::isLepOnMDC(kine,(HMatrixCategory*)getGeantMDCCat())     )
	{
	    if (kine->getTotalMomentum()>=50.){
	    //	    HRichUtilFunc::dumpKine(kine);
		n++;
		dalitzleps.Add(kine); // store lepton for further analysis 
		histoLepton(kine,hmomthetaallsingle);
	    }
	}
	
    }
    nbLeptons+=n;
    return n;//number of found and stored dalitz leptons in event
    
    
}

void HGeantKineAna::histoLepton(HGeantKine* lep,TH2D* h)
{
    if (lep)
    {
  	Float_t ptot=lep->getTotalMomentum();
	Float_t theta,phi;
	HRichUtilFunc::calcParticleAnglesKine(lep,theta,phi);
	h->Fill(ptot,theta);

    }

}
void HGeantKineAna::histoLeptonPair(TObjArray& p,TH2D* h)
{
	histoLepton((HGeantKine*)p[0],h);
	histoLepton((HGeantKine*)p[1],h);

}
void HGeantKineAna::checkAcceptedDalitzLeptons(TObjArray& dalitzleps)
{
    // !!! TObjArrays are auto-expanding !!!

    Int_t nbfoundRings=0;
    Int_t nbfoundGEANTDalitzPairs=0;
    Int_t nbfoundDetTracks = 0;
    Int_t nbfoundDetPairs = 0;
    //Int_t nbfoundGEANTLeptons = dalitzleps.GetLast()+1; // unused
    Int_t nbfoundGEANTSingleDalitzLeptons = 0;
    Int_t nbMatchedPairs = 0;
    Int_t nbMatchedTracks = 0;
    Int_t nbMatchedGPairs = 0;
    Int_t nbMatchedGTracks = 0;
    // unlikely to find more than one pair in event
    // starting from zero, two consecutive leptons are always a pair
    TObjArray dalitzpairs(1);
    
    if (dalitzleps.GetLast()==0) // there is only one lepton !
    {
	// check lepton: has ring(s) ? has track(s) ?
	nbfoundGEANTSingleDalitzLeptons++;
	nbDalitzSingle++;
    }
    else
    {
	// search for pair
	nbfoundGEANTDalitzPairs = fillDalitzPairs(dalitzleps,dalitzpairs);
	nbDalitzPairs+=nbfoundGEANTDalitzPairs;
	if (nbfoundGEANTDalitzPairs)
	{
	    histoLeptonPair(dalitzpairs,hmomthetaallpair);
	    // pair found

	}
	else {} // still no pair found
    }

    // array that holds for each lepton the corresponding 
    // ring whose recognized center is closer to the 
    // center of all reflected photons than for any other
    // ring
    TObjArray mindistrings(1);
    for (Int_t i=0;i<(dalitzleps.GetLast()+1);i++)
    {
	HRichHitSim *r=checkRings((HGeantKine*)dalitzleps[i]);
	if (r)
	{
	    nbfoundRings++;
	    mindistrings.Add(r);
	}
    }
    nbRings+=nbfoundRings;
    if (nbfoundRings)
    {
	// find all tracks that were found in the detectors w/o 
	// additional information from GEANT and that end in
	// the respective ring
	// array to hold track arrays for each ring
	TObjArray foundTracks(1);
	nbfoundDetTracks = getTracksForRings(mindistrings,foundTracks);
	nbTracks+=nbfoundDetTracks;
	nbMatchedTracks = matchGEANTwithFoundTracks(dalitzleps,foundTracks);
	nbgMatchedTracks+=nbMatchedTracks;

	nbMatchedGTracks = matchFoundTrackswithGEANT(dalitzleps,foundTracks);
	nbgMatchedGTracks+=nbMatchedGTracks;
	if (nbfoundDetTracks)
	{
	    // find all pairs formed with the above found tracks in this event
	    TObjArray foundPairs(1);
	    nbfoundDetPairs = getPairsForTracks(foundTracks,foundPairs);
	    nbTrackPairs+=nbfoundDetPairs;
	    nbMatchedPairs = matchGEANTwithFoundPair(dalitzpairs,foundPairs);
	    nbgMatchedPairs+=nbMatchedPairs;
	    nbMatchedGPairs = matchFoundPairwithGEANT(dalitzpairs,foundPairs);
	    nbgMatchedGPairs+=nbMatchedGPairs;
	}// endif found det tracks
    }//end if found rings

    // print statistics about found objects to screen
    

//      cout<<"Evt rings:"<<nbfoundRings<<" leptons:"<<nbfoundGEANTLeptons<<" dal pairs:"<<nbfoundGEANTDalitzPairs
//  	<<" tracks:"<<nbfoundDetTracks<<" pairs:"<<nbfoundDetPairs
//  	<<" singles:"<<nbfoundGEANTSingleDalitzLeptons<<endl;
//      cout<<"Sum rings:"<<nbRings<<" leptons:"<<nbLeptons<<" dal pairs:"<<nbDalitzPairs
//  	<<" tracks:"<<nbTracks<<" pairs:"<<nbTrackPairs
//  	<<" singles:"<<nbDalitzSingle<<endl;
//      cout<<"Nb GEANT confirmed reconstructed  pairs -- in Evt:"<<nbMatchedPairs
//  	<<"  Sum:"<<nbgMatchedPairs<<endl;
//      cout<<"Nb GEANT confirmed reconstructed tracks -- in Evt:"<<nbMatchedTracks
//  	<<"  Sum:"<<nbgMatchedTracks<<endl;
//      cout<<"Nb reconstructed  GEANT pairs -- in Evt:"<<nbMatchedGPairs
//  	<<"  Sum:"<<nbgMatchedGPairs<<endl;
//      cout<<"Nb reconstructed GEANT leptons -- in Evt:"<<nbMatchedGTracks
//  	<<"  Sum:"<<nbgMatchedGTracks<<endl;
}

Int_t HGeantKineAna::matchGEANTwithFoundTracks(TObjArray& Gl,TObjArray& Rt)
{
    // Gl: GEANT leptons
    // Rt: reconstructed tracks


    Int_t confirmedtracks = 0;
    Int_t *trackinds=new Int_t[Rt.GetLast()+1];
    for (Int_t kk=0;kk<(Rt.GetLast()+1);kk++) trackinds[kk]=-2;
    if ((Gl.GetLast()+1)>=1)
    {
	for (Int_t i=0;i<(Gl.GetLast()+1);i++)
	{
	    Int_t aTrack,aID;
	    ((HGeantKine*)Gl[i])->getParticle(aTrack, aID);
	    if ((Rt.GetLast()+1)>=1)
	    {
		for (Int_t j=0;j<(Rt.GetLast()+1);j++)
		{
		    HTrackInfo *t = ((HHitMatchSim*)Rt[j])->getTrackInfoObj();
		    for (Int_t k=0;k<MAXPARTICLES;k++)
		    {
			// count each track only once, even if
			// it contains more than one lepton
			if (t->getTrkNb(k)==aTrack &&
			    HRichCut::isNewIndex(getHitMatchCat()->
						 getIndex((HHitMatchSim*)Rt[j]),
						 trackinds,Rt.GetLast()+1) 
			    )
			{
			    confirmedtracks++;
			    break;
			}
		    }
		}
	    }
	}
    }
    delete [] trackinds;
    return confirmedtracks;
}
Int_t HGeantKineAna::matchFoundTrackswithGEANT(TObjArray& Gl,TObjArray& Rt)
{
    // Gl: GEANT leptons
    // Rt: reconstructed tracks


    Int_t confirmedlepton = 0;
    Int_t *trackinds=new Int_t[Rt.GetLast()+1];
    for (Int_t kk=0;kk<(Rt.GetLast()+1);kk++) trackinds[kk]=-2;
    if ((Gl.GetLast()+1)>=1)
    {
	for (Int_t i=0;i<(Gl.GetLast()+1);i++)
	{
	    Int_t aTrack,aID;
	    ((HGeantKine*)Gl[i])->getParticle(aTrack, aID);
	    if ((Rt.GetLast()+1)>=1)
	    {
		Bool_t lepisSeen=kFALSE;
		for (Int_t j=0;j<(Rt.GetLast()+1);j++)
		{
		    HTrackInfo *t = ((HHitMatchSim*)Rt[j])->getTrackInfoObj();
		    for (Int_t k=0;k<MAXPARTICLES;k++)
		    {
			// accept lepton only if it has its own track
			// the same track for two diff leps is forbidden
			// as this doesn't allow to reconstr and separate
			// the particle
			if (!lepisSeen &&
			    t->getTrkNb(k)==aTrack &&
			    HRichCut::isNewIndex(getHitMatchCat()->
						 getIndex((HHitMatchSim*)Rt[j]),
						 trackinds,Rt.GetLast()+1) 

			    )
			{
			    lepisSeen=kTRUE;
			    confirmedlepton++;
			    //Sl.Add((HGeantKine*)Gl[i]); // store reconstr. lepton
			    histoLepton((HGeantKine*)Gl[i],hmomthetareconsingle);
			    break;
			}
		    }
		}
	    }
	}
    }
    delete [] trackinds;
    return confirmedlepton;
}
Int_t HGeantKineAna::matchGEANTwithFoundPair(TObjArray& Gp,TObjArray& Rp)
{
    // Gp: GEANT pairs
    // Rp: reconstructed pairs
    
    Int_t confirmedpair = 0;

	if ((Gp.GetLast()+1)==2) // simple case with one pair
	{
	    // retrieve track numbers of leptons in pair
	    Int_t aTrack1, aTrack2;
	    Int_t aID1, aID2;
	    ((HGeantKine*)Gp[0])->getParticle(aTrack1, aID1);
	    ((HGeantKine*)Gp[1])->getParticle(aTrack2, aID2);
	    // loop over found pairs
	    for (Int_t i=0;i<(Rp.GetLast()+1);i++)
	    {
		// get indexes of reconst. tracks in pair
		Int_t ind1=((HDiHitMatch*)Rp[i])->getIndTrk1();
		Int_t ind2=((HDiHitMatch*)Rp[i])->getIndTrk2();
		// get tracks
		HHitMatchSim *trk1 = (HHitMatchSim*)(getHitMatchCat()
						     ->getObject(ind1));
		HHitMatchSim *trk2 = (HHitMatchSim*)(getHitMatchCat()
						     ->getObject(ind2));

		// get simulation part of reconstructed track
		HTrackInfo *t1 = trk1->getTrackInfoObj();
		HTrackInfo *t2 = trk2->getTrackInfoObj();

		Bool_t found11=kFALSE;
		Bool_t found12=kFALSE;
		Bool_t found21=kFALSE;
		Bool_t found22=kFALSE;
		// compare up-to MAXPARTICLES trk numbers
		// per reconstructed track to the given two
		for (Int_t k=0;k<MAXPARTICLES;k++)
		{
		    // if both trk numbers are found in one
		    // and none in the other, it is a Double_t 		    // with a fake !! exclude this pair !!
		    if (t1->getTrkNb(k)==aTrack1) found11=kTRUE;
		    if (t1->getTrkNb(k)==aTrack2) found12=kTRUE;
		    if (t2->getTrkNb(k)==aTrack1) found21=kTRUE;
		    if (t2->getTrkNb(k)==aTrack2) found22=kTRUE;
		}
		// good combinations
		if (
		    !( (found11&&found12&&!found21&&!found22) ||
		       (found21&&found22&&!found11&&!found12)
		     )   
		    )   
		{
		    // pair is confirmed now !!
		    confirmedpair++;
		    // here the pair could be added to arr !!
		}
		else
		{
		    trk1->dumpToStdoutSim();
		    trk2->dumpToStdoutSim();
		}
		
	    }//endfor
	}
	else if((Gp.GetLast()+1)<2)
	{
	    //no GEANT pair found, possible recon. fake pair
	}
	else if((Gp.GetLast()+1)>2)
	{
	    //more complicated, for Dalitz not needed
	}
	

	return confirmedpair;


}
Int_t HGeantKineAna::matchFoundPairwithGEANT(TObjArray& Gp,TObjArray& Rp)
{
    // Gp: GEANT pairs
    // Rp: reconstructed pairs
    
    Int_t confirmedpair = 0;

	if ((Gp.GetLast()+1)==2) // simple case with one pair
	{
	    // retrieve track numbers of leptons in pair
	    Int_t aTrack1, aTrack2;
	    Int_t aID1, aID2;
	    ((HGeantKine*)Gp[0])->getParticle(aTrack1, aID1);
	    ((HGeantKine*)Gp[1])->getParticle(aTrack2, aID2);
	    // loop over reconstructed pairs
	    for (Int_t i=0;i<(Rp.GetLast()+1);i++)
	    {
		// get indexes of reconst. tracks in pair
		Int_t ind1=((HDiHitMatch*)Rp[i])->getIndTrk1();
		Int_t ind2=((HDiHitMatch*)Rp[i])->getIndTrk2();
		// get tracks
		HHitMatchSim *trk1 = (HHitMatchSim*)(getHitMatchCat()
						     ->getObject(ind1));
		HHitMatchSim *trk2 = (HHitMatchSim*)(getHitMatchCat()
						     ->getObject(ind2));

		// get simulation part of reconstructed track
		HTrackInfo *t1 = trk1->getTrackInfoObj();
		HTrackInfo *t2 = trk2->getTrackInfoObj();

		Bool_t found11=kFALSE;
		Bool_t found12=kFALSE;
		Bool_t found21=kFALSE;
		Bool_t found22=kFALSE;
		// compare up-to MAXPARTICLES trk numbers
		// per reconstructed track to the given two from the leps
		for (Int_t k=0;k<MAXPARTICLES;k++)
		{
		    // if both trk numbers are found in one
		    // and none in the other, it is a Double_t 		    // with a fake !! exclude this pair !!
		    if (t1->getTrkNb(k)==aTrack1) found11=kTRUE;
		    if (t1->getTrkNb(k)==aTrack2) found12=kTRUE;
		    if (t2->getTrkNb(k)==aTrack1) found21=kTRUE;
		    if (t2->getTrkNb(k)==aTrack2) found22=kTRUE;
		}
		// good combinations
		if (
		    !( (found11&&found12&&!found21&&!found22) ||
		       (found21&&found22&&!found11&&!found12)
		     )   
		    )   
		{
		    // pair is confirmed now !!
		    confirmedpair++;
		    histoLeptonPair(Gp,hmomthetareconpair);
		    break; // stop looping over reconstructed pairs
		    // as the "real" pair was already found
		    // here the pair could be added to an arr !!
		}
		else
		{
		    trk1->dumpToStdoutSim();
		    trk2->dumpToStdoutSim();
		}
		
	    }//endfor reconstr. pairs
	}
	else if((Gp.GetLast()+1)<2)
	{
	    //no GEANT pair found, possible recon. fake pair
	}
	else if((Gp.GetLast()+1)>2)
	{
	    //more complicated, for Dalitz not needed
	}
	

	return confirmedpair;


}

Int_t HGeantKineAna::getPairsForTracks(TObjArray& tracks,TObjArray& pairs)
{

    HDiHitMatch *h=0;
    pIterDiHitMatch->Reset();
    Int_t nbFoundPairs=0;
    while(( h= (HDiHitMatch *)pIterDiHitMatch->Next())) 
    {
	//h->dumpToStdout();
	// indexes that made pair
	Int_t ind1=h->getIndTrk1();
	Int_t ind2=h->getIndTrk2();
	Int_t nbfoundtracksforpair=0;
	Int_t *trkinds = new Int_t[tracks.GetLast()+1];
	for (Int_t kk=0;kk<(tracks.GetLast()+1);kk++)trkinds[kk]=-2; 
//  	cout<<"ind1: "<<ind1<<"  ind2: "<<ind2<<endl;
//  	cout<<"track index: ";
	for (Int_t i=0;i<(tracks.GetLast()+1);i++)
	{
	    //	    ((HHitMatchSim*)tracks[i])->dumpToStdout();
	    // track index in track category
	    Int_t trkind=getHitMatchCat()->getIndex((HHitMatchSim*)tracks[i]);
	    //cout<<trkind<<" , ";
	    // check only against an index which did not occur before 
	    
	    if(
	       (trkind==ind1 || trkind==ind2) &&
	       HRichCut::isNewIndex(trkind,trkinds,tracks.GetLast()+1)
	       ) 
	    {
		//cout<<"found: "<<trkind<<" "<<ind1<<" "<<ind2<<endl;
		nbfoundtracksforpair++;
	    
	    
		if (nbfoundtracksforpair==2) 
		{
		    //cout<<"pair found"<<endl;
		    nbFoundPairs++;
		    pairs.Add(h);
		    break;
		}
	    }
	}
	//cout<<"pair loop"<<endl;
	//cout<<endl;
	delete [] trkinds;

    }
    return nbFoundPairs;
}
Int_t HGeantKineAna::getTracksForRings(TObjArray& rings,TObjArray& tracks)
{
    Int_t n=0;
    for (Int_t i=0;i<(rings.GetLast()+1);i++)
    {
	HHitMatchSim *track =0;
	pIterMatchHit->Reset();
	Int_t richind = getRichHitCat()->getIndex((HRichHitSim*)rings[i]);

	while((track=(HHitMatchSim *)pIterMatchHit->Next())!=0)
	{
	    if (
		HRichCut::isGoodRing((HRichHitSim*)rings[i]) &&
		HRichCut::isGoodTrack(track) &&
		track->getRichInd()==richind
		) 
	    {
		n++;
		tracks.Add(track);
	    }
	    
	}
    }
    return n;
}
Int_t HGeantKineAna::fillDalitzPairs(TObjArray& dalitzleps,TObjArray& dalitzpairs)
{
    Int_t n=0;
    for (Int_t i=0;i<(dalitzleps.GetLast()+1);i++)
    {
	for (Int_t j=i+1;j<(dalitzleps.GetLast()+1);j++) 
	{
	    Int_t aPar1, aMed1, aMech1;
	    Int_t aPar2, aMed2, aMech2;
	    Int_t aTrack1, aTrack2;
	    Int_t aID1, aID2;
	    
	    ((HGeantKine*)dalitzleps[i])->getParticle(aTrack1, aID1);
	    ((HGeantKine*)dalitzleps[j])->getParticle(aTrack2, aID2);
	    
	    ((HGeantKine*)dalitzleps[i])->getCreator(aPar1, aMed1, aMech1);
	    ((HGeantKine*)dalitzleps[j])->getCreator(aPar2, aMed2, aMech2);
	    
	    if (aPar1==aPar2 && aMed1==aMed2 && aMech1==aMech2)
	    {
		// pair found
		dalitzpairs.Add((HGeantKine*)dalitzleps[i]);
		dalitzpairs.Add((HGeantKine*)dalitzleps[j]);
		n++; //count pair
	    }
	}
    }
    return n;
}
Int_t HGeantKineAna::checkTrackPairs(TObjArray &t, TObjArray &p)
{
     
    HDiHitMatch *h=0;
    pIterDiHitMatch->Reset();
    Int_t nbFoundPairs=0;
    while(( h= (HDiHitMatch *)pIterDiHitMatch->Next())) 
    {
	Int_t ind1=h->getIndTrk1();
	Int_t ind2=h->getIndTrk2();
	Int_t foundtracks=0;
	Int_t *trkinds = new Int_t[t.GetLast()+1];
	for (Int_t j=0;j<(t.GetLast()+1);j++) trkinds[j]=-2;
	for (Int_t i=0;i<(t.GetLast()+1);i++)
	{
	    Int_t trkind=getHitMatchCat()->getIndex(t[i]);
	    if ((trkind==ind1 || trkind==ind2)&&
		HRichCut::isNewIndex(trkind,trkinds,t.GetLast()+1)) foundtracks++;

	}
	delete [] trkinds;
	if (foundtracks>=2) {nbFoundPairs++;p.Add(h);}
    }
    return nbFoundPairs;
}
Int_t HGeantKineAna::checkRingCorrelation(HRichHitSim* r,TObjArray& tracks)
{
    
    HHitMatchSim *track =0;
    pIterMatchHit->Reset();
    Int_t richind = getRichHitCat()->getIndex(r);
    Int_t n=0;
    while((track=(HHitMatchSim *)pIterMatchHit->Next())!=0)
    {
	if (
	    HRichUtilFunc::isGoodRing(r) &&
	    HRichUtilFunc::isGoodTrack(track) &&
	    track->getRichInd()==richind
	    ) 
	{
	    n++;
	    tracks.Add(track);
	}
	
    }
    return n;
}


Int_t HGeantKineAna::checkForPair(TObjArray& l,const Char_t *swt)
{
    if (l.GetLast()==0)
    {
	cout<<" **** pi0-Dalitz single lepton **** "<<endl;
	HRichUtilFunc::dumpKine((HGeantKine*)l[0]);
    }
    for (Int_t i=0;i<(l.GetLast()+1);i++)
    {
	for (Int_t j=i+1;j<(l.GetLast()+1);j++)
	{
	
	    Int_t aPar1, aMed1, aMech1;
	    Int_t aPar2, aMed2, aMech2;
	    Int_t aTrack1, aTrack2;
	    Int_t aID1, aID2;
	    
	    ((HGeantKine*)l[i])->getParticle(aTrack1, aID1);
	    ((HGeantKine*)l[j])->getParticle(aTrack2, aID2);

	    ((HGeantKine*)l[i])->getCreator(aPar1, aMed1, aMech1);
	    ((HGeantKine*)l[j])->getCreator(aPar2, aMed2, aMech2);

	    if (aPar1==aPar2 && aMed1==aMed2 && aMech1==aMech2)
	    {
		Float_t opang = HRichUtilFunc::calcOpeningAngleKine((HGeantKine*)l[i],(HGeantKine*)l[j]);
		Float_t invMass=HRichUtilFunc::calcInvMassKine((HGeantKine*)l[i],(HGeantKine*)l[j]);
		HGeantKine *gamma = 0;
		HGeantKine *pi0 = 0;
		HGeantKine *gamma1 =0;
		TString mode(swt);
		if (!mode.CompareTo("dalitz"))
		{
		// dalitz case
		    gamma = HRichUtilFunc::
			getPionDalitzGamma((HGeantKine*)l[i],
					   (HLinearCategory*)getGeantKineCat());
		    
		    if (gamma) pi0 =  HRichUtilFunc::
			findParent(gamma,(HLinearCategory*)getGeantKineCat());
		}
		else if (!mode.CompareTo("conv"))
		{
		// conversion case
		    gamma1 =  HRichUtilFunc::
			findParent((HGeantKine*)l[i],
				   (HLinearCategory*)getGeantKineCat());
		    if (gamma1) gamma = HRichUtilFunc::
			getSecondPionDecayGamma(gamma1,
						(HLinearCategory*)getGeantKineCat());
		     if (gamma1) pi0 =  HRichUtilFunc::
			findParent(gamma1,(HLinearCategory*)getGeantKineCat());
		}

		Bool_t validDecay = kFALSE;
		if (!mode.CompareTo("conv") && gamma1) validDecay=kTRUE;
		else if (!mode.CompareTo("dalitz")) validDecay = kTRUE;


		if (gamma&&pi0&&validDecay){
		
		// pi0 momentum
		Float_t ppi0X,ppi0Y,ppi0Z;
		pi0->getMomentum(ppi0X,ppi0Y,ppi0Z);


		// vector sum of leptons
		Float_t plep1X,plep1Y,plep1Z;
		Float_t plep2X,plep2Y,plep2Z;
		((HGeantKine*)l[i])->getMomentum(plep1X,plep1Y,plep1Z);
		((HGeantKine*)l[j])->getMomentum(plep2X,plep2Y,plep2Z);
		Float_t plepX = plep1X+plep2X;
		Float_t plepY = plep1Y+plep2Y;
		Float_t plepZ = plep1Z+plep2Z;
		Float_t pleptotal = TMath::Sqrt( plepX*plepX+
						 plepY*plepY+
						 plepZ*plepZ );
		HGeomVector psumleptons;
		psumleptons.setX(plepX);
		psumleptons.setY(plepY);
		psumleptons.setZ(plepZ);
		Float_t totMom1 = ((HGeantKine*)l[i])->getTotalMomentum();
		Float_t totMom2 = ((HGeantKine*)l[j])->getTotalMomentum();
		hmompairkine->Fill(totMom1,totMom2);
		//		Float_t momgamma1 = gamma1->getTotalMomentum();


		Float_t theta1,phi1;
		Float_t theta2,phi2;
		HRichUtilFunc::calcParticleAnglesKine((HGeantKine*)l[i],theta1,phi1);
		HRichUtilFunc::calcParticleAnglesKine((HGeantKine*)l[j],theta2,phi2);
		hthetapairkine->Fill(theta1,theta2);
		// take care of gamma
		Float_t momgamma = gamma->getTotalMomentum();
		HGeomVector pgamma;
		Float_t pgammaX,pgammaY,pgammaZ;
		
		gamma->getMomentum(pgammaX,pgammaY,pgammaZ);
		pgamma.setX(pgammaX);
		pgamma.setY(pgammaY);
		pgamma.setZ(pgammaZ);
		
		// take care of second gamma
		//  if (!mode.CompareTo("conv"))
//  		{
//  		    Float_t momgamma1 = gamma1->getTotalMomentum();
//  		    Float_t theta_gamma,phi_gamma;
//  		    Float_t theta_gamma1,phi_gamma1;
		    
//  		    HRichUtilFunc::calcParticleAnglesKine(gamma,theta_gamma,phi_gamma);
//  		    HRichUtilFunc::calcParticleAnglesKine(gamma1,theta_gamma1,phi_gamma1);
//  		    hgammathetatheta->Fill(theta_gamma,theta_gamma1);
//  		    hgammamommom->Fill(momgamma,momgamma1);
//  		}

		// check momentum conservation
		Float_t psumX,psumY,psumZ;
		psumX = plep1X+plep2X+pgammaX-ppi0X;
		psumY = plep1Y+plep2Y+pgammaY-ppi0Y;
		psumZ = plep1Z+plep2Z+pgammaZ-ppi0Z;
		
		// calc opening angle between gamma mom 3vector 
		// and sum of lepton mom 3vectors 
		Float_t opanglepsgamma = HRichUtilFunc::calcOpeningAngleVectors(pgamma,psumleptons);
		// calc inv mass of decayed pion0 (134.976 MeV, c.f.:data booklet)
		Float_t invMassPion = 
		    2.*
		    sin( (1./HConstant::rad2deg()) *opanglepsgamma/2. )*
		    sqrt(momgamma*pleptotal);


		// dump to screen
//  		if (!mode.CompareTo("dalitz"))
//  		    cout<<" **** pi0-Dalitz decay (";
//  		else if (!mode.CompareTo("conv"))
//  		    cout<<" **** pi0->2-gamma decay (";
//  		cout<<"opang:"<<opang
//  		    <<"; lep pair inv. mass:"<<invMass;
//  		//cout<<"; trk nb1:"<<aTrack1<<"; trk nb2:"<<aTrack2;
//  		cout<<")"<<endl;
//    		cout<<" **** mom_gamma:"<<momgamma
//    		    <<" pion mass:"<<invMassPion<<endl;

  		//cout<<" **** mom sum 3-vector p("<<psumX<<","<<psumY<<","<<psumZ<<")  "

//              HRichUtilFunc::dumpKine((HGeantKine*)l[i]);
//  		HRichUtilFunc::dumpKine((HGeantKine*)l[j]);
//  		HRichUtilFunc::dumpKine(gamma);
//  		HRichUtilFunc::dumpKine(pi0);
		hpi0invmass->Fill(invMassPion);
		hdalitzopang->Fill(opang);
		hdalitzinvmass->Fill(invMass);
		}
	    }
	}

    }
    return 1;
}

void HGeantKineAna::calcParticleAnglesV(HGeantKine *kine,Float_t &theta, Float_t &phi)
{
    // input kine object with momentum vector
    // output theta and phi of trajectory
    Float_t xMom,yMom,zMom;
    kine->getMomentum(xMom,yMom,zMom);
    HGeomVector2 vec;
    vec.setX(xMom);
    vec.setY(yMom);
    vec.setZ(zMom);
    
    vec/=vec.length();//norm

    Float_t rad;
    vec.sphereCoord(rad,theta,phi);
    //cout<<"theta: "<<theta<<" phi: "<<phi<<endl;
}

Float_t HGeantKineAna::calcOpeningAngleV(HGeantKine *kine1,HGeantKine *kine2)
{
    //input 2 kine objects
    //output opening angle of trajectories
    Float_t rad2deg = 180./TMath::Pi();
    HGeomVector vec1;
    if (kine1){
    Float_t xMom1,yMom1,zMom1;
    kine1->getMomentum(xMom1,yMom1,zMom1);
    vec1.setX(xMom1);
    vec1.setY(yMom1);
    vec1.setZ(zMom1); 
    
    vec1/=vec1.length();//norm
    }
    HGeomVector vec2;
    if (kine2){
    Float_t xMom2,yMom2,zMom2;
    kine2->getMomentum(xMom2,yMom2,zMom2);
    vec2.setX(xMom2);
    vec2.setY(yMom2);
    vec2.setZ(zMom2);
    
    vec2/=vec2.length();//norm
    }
    
    Float_t cosfOpeningAngle = vec1.scalarProduct(vec2);
    //cout<<cosfOpeningAngle<<endl;
    Float_t fOpeningAngle=0.;
    if (TMath::Abs(cosfOpeningAngle) <= 1) 
	fOpeningAngle=TMath::ACos(cosfOpeningAngle) * rad2deg;
    
    return fOpeningAngle;
}

void HGeantKineAna::calcParticleAnglesT(HGeantKine *kine,Float_t &fpt, Float_t &fpp)
{
    //input kine object
    //output theta and phi of trajectory

    //dumpKine(kine);
    Float_t rad2deg = 180./TMath::Pi();
    Float_t pX,pY,pZ;
    pX=pY=pZ=0;
    kine->getMomentum(pX,pY,pZ);
    Float_t pT = TMath::Sqrt(pX*pX+pY*pY);
    Float_t pTot = kine->getTotalMomentum();
    fpt = TMath::ASin(pT/pTot) * rad2deg; //theta (polar) of particle
    fpp = TMath::ASin(pY/pT) *rad2deg; //phi (azimuth) of particle
    if (pX<0.) fpp = 180.-fpp;
    else if (pY<0.) fpp+=360.;
    //cout<<" Geant theta: "<<fpt
    //<< "Geant phi: "<<fpp<<endl;

}

Float_t HGeantKineAna::calcOpeningAngleT(Float_t t1,Float_t p1,
				      Float_t t2,Float_t p2)
{
    //input theta and phi of two trajectories
    //returns opening angle 
    Float_t rad2deg = 180./TMath::Pi();
    HGeomVector vec1,vec2;

    vec1.setX(TMath::Sin(t1/rad2deg) * TMath::Cos(p1/rad2deg));
    vec1.setY(TMath::Sin(t1/rad2deg) * TMath::Sin(p1/rad2deg));
    vec1.setZ(TMath::Cos(t1/rad2deg));

    vec2.setX(TMath::Sin(t2/rad2deg) * TMath::Cos(p2/rad2deg));
    vec2.setY(TMath::Sin(t2/rad2deg) * TMath::Sin(p2/rad2deg));
    vec2.setZ(TMath::Cos(t2/rad2deg));

    
    Float_t cosfOpeningAngle = vec1.scalarProduct(vec2);
    cout<<cosfOpeningAngle<<endl;
    Float_t fOpeningAngle=0.;
    if (TMath::Abs(cosfOpeningAngle) <= 1) 
	fOpeningAngle=TMath::ACos(cosfOpeningAngle) * rad2deg;
    
    return fOpeningAngle;
}

void HGeantKineAna::dumpKine(HGeantKine *kine)
{
    if (kine)
    {
    Int_t aTrack, aID;
    Int_t aPar, aMed, aMech;
    Float_t ax, ay, az;
    Float_t apx, apy, apz;
    Float_t aInfo, aWeight;
    Float_t ptot;
    kine->getParticle(aTrack,aID);
    kine->getCreator(aPar,aMed,aMech);
    kine->getVertex(ax,ay,az);
    kine->getMomentum(apx,apy,apz);
    kine->getGenerator(aInfo,aWeight);
    ptot=kine->getTotalMomentum();
    cout<<"-----------------------------------------------------------"<<endl;
    cout<<"--- GEANT track number: "<<aTrack<<endl;
    cout<<"--- track number of parent particle: "<<aPar<<endl;
    cout<<"--- GEANT particle ID: "<<aID<<endl;
    cout<<"--- GEANT medium of creation: "<<aMed<<endl;
    cout<<"--- GEANT creation mechanism: "<<aMech<<endl;
    cout<<"--- x of track vertex (in mm): "<<ax<<endl;
    cout<<"--- y                        : "<<ay<<endl;
    cout<<"--- z                        : "<<az<<endl;
    cout<<"--- x component of particle momentum (in MeV/c): "<<apx<<endl;
    cout<<"--- y                                          : "<<apy<<endl;
    cout<<"--- z                                          : "<<apz<<endl;
    cout<<"--- total momentum                             : "<<ptot<<endl;
    cout<<"--- event generator info: "<<aInfo<<endl;
    cout<<"--- associated weight: "<<aWeight<<endl;
    cout<<"--- "<<endl;
    cout<<"--- "<<endl;
    cout<<"--- "<<endl;
    cout<<"--- "<<endl;
    cout<<"-----------------------------------------------------------"<<endl;
    }
}
void HGeantKineAna::histoKine(HGeantKine *kine)
{
    if (kine)
    {
	Int_t aTrack, aID;
	Int_t aPar, aMed, aMech;
	Float_t ax, ay, az;
	Float_t apx, apy, apz;
	Float_t aInfo, aWeight;
	Float_t ptot;
	kine->getParticle(aTrack,aID);
	kine->getCreator(aPar,aMed,aMech);
	kine->getVertex(ax,ay,az);
	kine->getMomentum(apx,apy,apz);
	kine->getGenerator(aInfo,aWeight);
  	ptot=kine->getTotalMomentum();
	Float_t theta,phi;
	HRichUtilFunc::calcParticleAnglesKine(kine,theta,phi);
	hidkine->Fill(aID);
	hmedkine->Fill(aMed);
	hmechkine->Fill(aMech);
	hmomkine->Fill(ptot);
	hmomthetakine->Fill(theta,ptot);
	hparentkine->Fill(aPar);
	hvertexkine->Fill(ax,ay,az);
    }
}

HGeantKine* HGeantKineAna::findParent(HGeantKine *kine)
{
    Int_t aPar, aMed, aMech;
    kine->getCreator(aPar,aMed,aMech);
    if (aPar){
	HIterator* iter_kine2 = (HIterator*)(((HLinearCategory*)
					      getGeantKineCat())
					     ->MakeIterator("native"));
	iter_kine2->Reset();
	HGeantKine *kine2=0;
	Int_t aTrackParent,aIDParent;
	while((kine2=(HGeantKine *)iter_kine2->Next())!=0)
	{
	    kine2->getParticle(aTrackParent,aIDParent);;
	    if (aPar == aTrackParent) 
	    {
		//if (kDumpIt) dumpKine(findParent(kine2));//recursive research for relatives
		return kine2;
	    }
	}
    }
    return 0;
}

Bool_t HGeantKineAna::finalize() {
    //cout<<"in finalize"<<endl;
//      TFile ff("largedistevtlist.root","RECREATE");
//      ff.cd();
//      el->Print("all");
//      el->Write();
//      ff.Close();
    HRichUtilFunc::saveHistos(pFileOut,pHistArray);
    pFileOut->Close();

    cout<<"Sum rings:"<<nbRings<<" leptons:"<<nbLeptons<<" dal pairs:"<<nbDalitzPairs
	<<" tracks:"<<nbTracks<<" pairs:"<<nbTrackPairs
	<<" singles:"<<nbDalitzSingle<<endl;
    cout<<"Nb GEANT confirmed reconstructed pairs -- Sum:"<<nbgMatchedPairs<<endl;
    cout<<"Nb of reconstructed GEANT pairs -- Sum:"<<nbgMatchedGPairs<<endl;
    cout<<"Nb GEANT confirmed reconstructed tracks -- Sum:"<<nbgMatchedTracks<<endl;
    cout<<"Nb of reconstructed GEANT leptons -- Sum:"<<nbgMatchedGTracks<<endl;
    //    delete el;
    return kTRUE;
}

void HGeantKineAna::saveHistos()
{
    pFileOut->cd();
    // write histograms
    for (Int_t i=0;i<(pHistArray->GetLast()+1);i++)
    {
	( (TH1*)(*pHistArray)[i] )->Write();
    }
    //pHistElePosMomConv->Write();
    //pHistOpenAngle->Write();

}