// -------------------------------------------------------------------------
// -----                 PndMCTestLambdaLambdaBar source file             -----
// -----                  Created 18/07/08  by S.Esch              -----
// -------------------------------------------------------------------------
// libc includes
#include <iostream>

// Root includes
#include "TROOT.h"
#include "TClonesArray.h"
#include "TVector3.h"
#include "TLorentzVector.h"

// framework includes
#include "FairRootManager.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "FairHit.h"
#include "FairMultiLinkedData.h"
#include "FairEventHeader.h"
#include "PndMCTestLambdaLambdaBar.h"

#include "PndSdsMCPoint.h"
#include "PndSdsHit.h"
#include "PndGemHit.h"
#include "PndSttHit.h"
#include "PndMCEntry.h"
#include "PndMCTrack.h"
#include "PndGemMCPoint.h"
#include "PndTrackCand.h"

#include "PndDetectorList.h"
#include <iomanip>

// -----   Default constructor   -------------------------------------------
PndMCTestLambdaLambdaBar::PndMCTestLambdaLambdaBar()
	: FairTask("Creates PndMC test")
{
}
// -------------------------------------------------------------------------

// -----   Destructor   ----------------------------------------------------
PndMCTestLambdaLambdaBar::~PndMCTestLambdaLambdaBar()
{
}

// -----   Public method Init   --------------------------------------------
InitStatus PndMCTestLambdaLambdaBar::Init()
{
	eventcounter=0;
	fLogger->GetLogger();
	FairRootManager* ioman = FairRootManager::Instance();
  	if (!ioman) {
  		std::cout << "-E- PndMCTestLambdaLambdaBar::Init: "
  				<< "RootManager not instantiated!" << std::endl;
  		return kFATAL;
  	}

  	fMCMatch 		= 	(PndMCMatch*)ioman->GetObject("MCMatch");
  	fEventHeader 	= 	(TClonesArray*)ioman->GetObject("EventHeader.");
  	fMCPoint 		= 	(TClonesArray*)ioman->GetObject("MVDPoint");
  	fMCTrack 		= 	(TClonesArray*)ioman->GetObject("MCTrack");

  	fPositionDecayVertexLambdaLambdaBar2D = new TH2F("fPositionDecayVertexLambdaLambdaBar2D", "Position Decay Vertex #Lambda #bar{#Lambda} ; z Position [cm] ; R [cm]", 900,0,90,200,0,20);
    fPositionDecayVertexLambdaBar2D = new TH2F("fPositionDecayVertexLambdaBar2D", "Position Decay Vertex #bar{#Lambda} ; z Position [cm] ; R [cm]", 900,0,90,200,0,20);
    fPositionDecayVertexLambda2D = new TH2F("fPositionDecayVertexLambda2D", "Position Decay Vertex #Lambda ; z Position [cm] ; R [cm]", 900,0,90,200,0,20);
    fPositionBothLambdasDecayed = new TH2F("fPositionBothLambdasDecayed", "Position both #Lambda decayed; z Position [cm] ; R [cm]", 900,0,90,200,0,20);
    fMeanLifeLambdaCMS = new TH1F("fMeanLifeLambdaCMS", "Mean Life #Lambda CMS; #frac{d m_{0}}{p} [cm] ; # Events",400,0, 100);
    fMeanLifeLambdaBarCMS = new TH1F("fMeanLifeLambdaBarCMS", "Mean Life #bar{#Lambda} CMS; #frac{d m_{0}}{p} [cm] ; # Events", 400,0, 100);
    fMeanLifeLambdaLAB = new TH1F("fMeanLifeLambdaLAB", "Mean Life #Lambda LAB; d [cm] ; # Events", 400,0, 100);
    fMeanLifeLambdaBarLAB = new TH1F("fMeanLifeLambdaBarLAB", "Mean Life #bar{#Lambda} LAB; d [cm] ; # Events", 400,0, 100);
    fArmenterosPodolanskiPlot = new TH2F("fArmenterosPodolanskiPlot","Armenteros-Podolanski-Plot; #alpha; P_t [GeV/c]",400,-1.5,1.5,400,0,0.12);

    fAngularDistributionLambdaCMSCosTheta = new TH1F("fAngularDistributionLambdaCMSCosTheta", "Angular Distribution All Particles; Cos(#Theta); # ", 40,-1,1);
    fAngularDistributionLambdabarCMSCosTheta = new TH1F("fAngularDistributionLambdabarCMSCosTheta", "Angular Distribution All Particles; Cos(#Theta); # ", 40,-1,1);
    fAngularDistributionProtonCMSCosTheta = new TH1F("fAngularDistributionProtonCMSCosTheta", "Angular Distribution Proton ; Cos(#Theta); #", 40,-1,1);
    fAngularDistributionAntiprotonCMSCosTheta = new TH1F("fAngularDistributionAntiprotonCMSCosTheta", "Angular Distribution Antiproton; Cos(#Theta); #", 40,-1,1);
    fAngularDistributionPionplusCMSCosTheta = new TH1F("fAngularDistributionPionplusCMSCosTheta", "Angular Distribution Pionplus; Cos(#Theta); #", 40,-1,1);
    fAngularDistributionPionminusCMSCosTheta = new TH1F("fAngularDistributionPionminusCMSCosTheta", "Angular Distribution pionminus; Cos(#Theta); #", 40,-1,1);

    fAngularDistributionLambdaCMSTheta = new TH1F("fAngularDistributionLambdaCMSTheta", "Angular Distribution All Particles; #Theta; # ", 320,0,180);
    fAngularDistributionLambdabarCMSTheta = new TH1F("fAngularDistributionLambdabarCMSTheta", "Angular Distribution All Particles; #Theta; # ", 320,0,180);
    fAngularDistributionProtonCMSTheta = new TH1F("fAngularDistributionProtonCMSTheta", "Angular Distribution Proton ; #Theta; #", 320,0,180);
    fAngularDistributionAntiprotonCMSTheta = new TH1F("fAngularDistributionAntiprotonCMSTheta", "Angular Distribution Antiproton; #Theta; #", 320,0,180);
    fAngularDistributionPionplusCMSTheta = new TH1F("fAngularDistributionPionplusCMSTheta", "Angular Distribution Pionplus; #Theta; #",320,0,180);
    fAngularDistributionPionminusCMSTheta = new TH1F("fAngularDistributionPionminusCMSTheta", "Angular Distribution pionminus; #Theta; #", 320,0,180);

    fAngularDistributionLambdaLABCosTheta = new TH1F("fAngularDistributionLambdaLABCosTheta", "Angular Distribution All Particles; Cos(#Theta); # ", 40,-1,1);
    fAngularDistributionLambdabarLABCosTheta = new TH1F("fAngularDistributionLambdabarLABCosTheta", "Angular Distribution All Particles; Cos(#Theta); # ", 40,-1,1);
    fAngularDistributionProtonLABCosTheta = new TH1F("fAngularDistributionProtonLABCosTheta", "Angular Distribution Proton ; Cos(#Theta); #", 40,-1,1);
    fAngularDistributionAntiprotonLABCosTheta = new TH1F("fAngularDistributionAntiprotonLABCosTheta", "Angular Distribution Antiproton; Cos(#Theta); #", 40,-1,1);
    fAngularDistributionPionplusLABCosTheta = new TH1F("fAngularDistributionPionplusLABCosTheta", "Angular Distribution Pionplus; Cos(#Theta); #", 40,-1,1);
    fAngularDistributionPionminusLABCosTheta = new TH1F("fAngularDistributionPionminusLABCosTheta", "Angular Distribution pionminus; Cos(#Theta); #", 40,-1,1);

    fAngularDistributionLambdaLABTheta = new TH1F("fAngularDistributionLambdaLABTheta", "Angular Distribution All Particles; #Theta; # ", 320,0,180);
    fAngularDistributionLambdabarLABTheta = new TH1F("fAngularDistributionLambdabarLABTheta", "Angular Distribution All Particles; #Theta; # ", 320,0,180);
    fAngularDistributionProtonLABTheta = new TH1F("fAngularDistributionProtonLABTheta", "Angular Distribution Proton ; #Theta; #", 320,0,180);
    fAngularDistributionAntiprotonLABTheta = new TH1F("fAngularDistributionAntiprotonLABTheta", "Angular Distribution Antiproton; #Theta; #", 320,0,180);
    fAngularDistributionPionplusLABTheta = new TH1F("fAngularDistributionPionplusLABTheta", "Angular Distribution Pionplus; #Theta; #",320,0,180);
    fAngularDistributionPionminusLABTheta = new TH1F("fAngularDistributionPionminusLABTheta", "Angular Distribution pionminus; #Theta; #",320,0,180);

    fPzVsPtLambda = new TH2F("fPzVsPtLambda", "Pz vs. Pt #Lambda; Pz [GeV/c] ; Pt [GeV/c]", 200, 0, 2, 40, 0, 0.4);
    fPzVsPtLambdabar = new TH2F("fPzVsPtLambdabar", "Pz vs. Pt #bar{#Lambda} ; Pz [GeV/c] ; Pt [GeV/c]", 200, 0, 2, 40, 0, 0.4);

    fPositionMvdMcPoint = new TH2F("fPositionMvdMcPoint","Position MVD MCPoint;z [cm]; R [cm]",18000,0,90,4000,0,20);

    fDetectorCoverageMVDMCPointsProton = new TH2F("fDetectorCoverageMVDMCPointsProton", "Detector Coverage MVD Points Proton; #Phi ; #Theta", 360,-180,180,180,0,180);
    fDetectorCoverageMVDMCPointsAntiproton = new TH2F("fDetectorCoverageMVDMCPointsAntiproton", "Detector Coverage MVD Points Antiproton; #Phi ; #Theta", 360,-180,180,180,0,180);
    fDetectorCoverageMVDMCPointsPionplus= new TH2F("fDetectorCoverageMVDMCPointsPionplus", "Detector Coverage MVD Points Pionplus; #Phi ; #Theta", 360,-180,180,180,0,180);
    fDetectorCoverageMVDMCPointsPionminus= new TH2F("fDetectorCoverageMVDMCPointsPionminus", "Detector Coverage MVD Points Pionminus; #Phi ; #Theta", 360,-180,180,180,0,180);

    fMCTracksPerEvent = new TH1I("fMCTracksPerEvent","fMCTracksPerEvent;#FinalState Tracks/Event;#Events",7,-1.5,5.5);

    std::cout << "-I- PndMCTestLambdaLambdaBar::Init: Initialization successfull" << std::endl;

  return kSUCCESS;
}

// -------------------------------------------------------------------------
void PndMCTestLambdaLambdaBar::SetParContainers()
{
  // Get Base Container
//  FairRun* ana = FairRun::Instance();
//  FairRuntimeDb* rtdb=ana->GetRuntimeDb();
	  fGeoH = PndGeoHandling::Instance();
	  fGeoH->SetParContainers();
}

// -----   Public method Exec   --------------------------------------------
void PndMCTestLambdaLambdaBar::Exec(Option_t* opt)
{
	if(fVerbose>0)
	{
		std::cout << "============= Begin PndMCTestLambdaLambdaBar::Exec for Event " << eventcounter << std::endl;
	}

	FairRootManager* ioman = FairRootManager::Instance();

TVector3 momentumproton		;
TVector3 momentumantiproton	;
TVector3 momentumpionplus	;
TVector3 momentumpionminus	;

double positionproton 		=0;
double positionprotonx 		=0;
double positionprotony 		=0;
double positionprotonz 		=0;

double positionantiproton 		=0;
double positionantiprotonx 		=0;
double positionantiprotony 		=0;
double positionantiprotonz 		=0;

double positionpionplus 		=0;
double positionpionplusx 		=0;
double positionpionplusy 		=0;
double positionpionplusz 		=0;

double positionpionminus 		=0;
double positionpionminusx 		=0;
double positionpionminusy 		=0;
double positionpionminusz 		=0;

double momentumlambda		=0;
double momentumlambdabar	=0;

double alphalambda		=0;
double alphalambdabar	=0;

double ptproton		=0;
double ptantiproton	=0;
double ptpionplus	=0;
double ptpionminus	=0;

int protontrack		=0;
int antiprotontrack	=0;
int pionminustrack	=0;
int pionplustrack	=0;

TLorentzVector lorentzvectorproton;
TLorentzVector lorentzvectorantiproton;
TLorentzVector lorentzvectorpionplus;
TLorentzVector lorentzvectorpionminus;
TLorentzVector lorentzvectorlambda;
TLorentzVector lorentzvectorlambdabar;
TLorentzVector lorentzvectorgesamt;

PndSdsMCPoint* myMCPoint;

for(int t=0;t<fMCPoint->GetEntries();t++)
{
	myMCPoint=(PndSdsMCPoint*)fMCPoint->At(t);
	if(myMCPoint->GetTrackID() <4)
	{
	fPositionMvdMcPoint->Fill(myMCPoint->GetZ(),TMath::Sqrt((myMCPoint->GetY()*myMCPoint->GetY())+(myMCPoint->GetX()*myMCPoint->GetX())));
	}
}

PndMCTrack* myTrackVertex;

for(int w=0;w<fMCTrack->GetEntries() && w<10 ;w++)
{
	myTrackVertex=(PndMCTrack*)fMCTrack->At(w);
	int STTPoint = myTrackVertex->GetNPoints(kSTT);
	int GEMPoint = myTrackVertex->GetNPoints(kGEM);
	int FTSPoint = myTrackVertex->GetNPoints(kFTS);
	int MVDPoint = myTrackVertex->GetNPoints(kMVD);

	if(myTrackVertex->GetPdgCode()==2212 && myTrackVertex->GetMotherID() == -1)
	{
		std::cout <<" Protontrack " << "MVD " << MVDPoint <<  " STT " << STTPoint << " GEM " << GEMPoint <<  " FTS " << FTSPoint << std::endl;

		if(STTPoint!=0 || GEMPoint !=0 || FTSPoint!=0 || MVDPoint !=0 )
		{
			protontrack=1;

		}
		momentumproton = myTrackVertex->GetMomentum();
		ptproton = myTrackVertex->GetPt();
		positionproton = (myTrackVertex->GetStartVertex()).Mag();
		positionprotonx= (myTrackVertex->GetStartVertex()).X();
		positionprotony= (myTrackVertex->GetStartVertex()).Y();
		positionprotonz= (myTrackVertex->GetStartVertex()).Z();
		lorentzvectorproton=myTrackVertex->Get4Momentum();
	}
	else if(myTrackVertex->GetPdgCode()==-2212 && myTrackVertex->GetMotherID() == -1)
	{
		std::cout <<" Antiprotontrack " << "MVD " << MVDPoint <<  " STT " << STTPoint << " GEM " << GEMPoint <<  " FTS " << FTSPoint << std::endl;

		if(STTPoint!=0 || GEMPoint !=0 || FTSPoint!=0 || MVDPoint !=0 )
		{
			antiprotontrack=1;

		}
		momentumantiproton = myTrackVertex->GetMomentum();
		ptantiproton = myTrackVertex->GetPt();
		positionantiproton = (myTrackVertex->GetStartVertex()).Mag();
		positionantiprotonx= (myTrackVertex->GetStartVertex()).X();
		positionantiprotony= (myTrackVertex->GetStartVertex()).Y();
		positionantiprotonz= (myTrackVertex->GetStartVertex()).Z();
		lorentzvectorantiproton=myTrackVertex->Get4Momentum();
	}
	else if(myTrackVertex->GetPdgCode()==211 && myTrackVertex->GetMotherID() == -1)
	{
		std::cout <<" Pionplus track" << "MVD " << MVDPoint <<  " STT " << STTPoint << " GEM " << GEMPoint <<  " FTS " << FTSPoint << std::endl;

		if(STTPoint!=0 || GEMPoint !=0 || FTSPoint!=0 || MVDPoint !=0 )
		{
			pionplustrack=1;

		}
		momentumpionplus = myTrackVertex->GetMomentum();
		ptpionplus = myTrackVertex->GetPt();
		positionpionplus = (myTrackVertex->GetStartVertex()).Mag();
		positionpionplusx= (myTrackVertex->GetStartVertex()).X();
		positionpionplusy= (myTrackVertex->GetStartVertex()).Y();
		positionpionplusz= (myTrackVertex->GetStartVertex()).Z();
		lorentzvectorpionplus=myTrackVertex->Get4Momentum();
	}
	else if(myTrackVertex->GetPdgCode()==-211 && myTrackVertex->GetMotherID() == -1)
	{
		std::cout <<" Pionminus track " << "MVD " << MVDPoint <<  " STT " << STTPoint << " GEM " << GEMPoint <<  " FTS " << FTSPoint << std::endl;

		if(STTPoint!=0 || GEMPoint !=0 || FTSPoint!=0 || MVDPoint !=0 )
		{
			pionminustrack=1;

		}
		momentumpionminus = myTrackVertex->GetMomentum();
		ptpionminus = myTrackVertex->GetPt();
		positionpionminus = (myTrackVertex->GetStartVertex()).Mag();
		positionpionminusx= (myTrackVertex->GetStartVertex()).X();
		positionpionminusy= (myTrackVertex->GetStartVertex()).Y();
		positionpionminusz= (myTrackVertex->GetStartVertex()).Z();
		lorentzvectorpionminus=myTrackVertex->Get4Momentum();
	}
}
lorentzvectorlambda = lorentzvectorproton+lorentzvectorpionminus;
lorentzvectorlambdabar = lorentzvectorpionplus+lorentzvectorantiproton;
lorentzvectorgesamt = lorentzvectorproton+lorentzvectorpionminus+lorentzvectorpionplus+lorentzvectorantiproton;
TVector3 betavector = lorentzvectorgesamt.BoostVector();

fPzVsPtLambda->Fill(lorentzvectorlambda.Pz(),lorentzvectorlambda.Pt());
fPzVsPtLambdabar->Fill(lorentzvectorlambdabar.Pz(),lorentzvectorlambdabar.Pt());

fMCTracksPerEvent->Fill(protontrack+antiprotontrack+pionplustrack+pionminustrack);
std::cout <<"protontrack "<< protontrack << std::endl;
std::cout <<"antiprotontrack "<< antiprotontrack << std::endl;
std::cout <<"pionplustrack "<< pionplustrack << std::endl;
std::cout <<"pionminustrack "<< pionminustrack << std::endl;

fAngularDistributionLambdaLABCosTheta->Fill(lorentzvectorlambda.CosTheta());
fAngularDistributionLambdabarLABCosTheta->Fill(lorentzvectorlambdabar.CosTheta());
fAngularDistributionProtonLABCosTheta->Fill(lorentzvectorproton.CosTheta());
fAngularDistributionAntiprotonLABCosTheta->Fill(lorentzvectorantiproton.CosTheta());
fAngularDistributionPionplusLABCosTheta->Fill(lorentzvectorpionplus.CosTheta());
fAngularDistributionPionminusLABCosTheta->Fill(lorentzvectorpionminus.CosTheta());

fAngularDistributionLambdaLABTheta->Fill(lorentzvectorlambda.Theta()*TMath::RadToDeg());
fAngularDistributionLambdabarLABTheta->Fill(lorentzvectorlambdabar.Theta()*TMath::RadToDeg());
fAngularDistributionProtonLABTheta->Fill(lorentzvectorproton.Theta()*TMath::RadToDeg());
fAngularDistributionAntiprotonLABTheta->Fill(lorentzvectorantiproton.Theta()*TMath::RadToDeg());
fAngularDistributionPionplusLABTheta->Fill(lorentzvectorpionplus.Theta()*TMath::RadToDeg());
fAngularDistributionPionminusLABTheta->Fill(lorentzvectorpionminus.Theta()*TMath::RadToDeg());

fDetectorCoverageMVDMCPointsProton->Fill(lorentzvectorproton.Phi()*TMath::RadToDeg(),lorentzvectorproton.Theta()*TMath::RadToDeg());
fDetectorCoverageMVDMCPointsAntiproton->Fill(lorentzvectorantiproton.Phi()*TMath::RadToDeg(),lorentzvectorantiproton.Theta()*TMath::RadToDeg());
fDetectorCoverageMVDMCPointsPionplus->Fill(lorentzvectorpionplus.Phi()*TMath::RadToDeg(),lorentzvectorpionplus.Theta()*TMath::RadToDeg());
fDetectorCoverageMVDMCPointsPionminus->Fill(lorentzvectorpionminus.Phi()*TMath::RadToDeg(),lorentzvectorpionminus.Theta()*TMath::RadToDeg());

lorentzvectorlambda.Boost(-betavector);
lorentzvectorlambdabar.Boost(-betavector);
lorentzvectorantiproton.Boost(-betavector);
lorentzvectorproton.Boost(-betavector);
lorentzvectorpionplus.Boost(-betavector);
lorentzvectorpionminus.Boost(-betavector);

fAngularDistributionLambdaCMSCosTheta->Fill(lorentzvectorlambda.CosTheta());
fAngularDistributionLambdabarCMSCosTheta->Fill(lorentzvectorlambdabar.CosTheta());
fAngularDistributionProtonCMSCosTheta->Fill(lorentzvectorproton.CosTheta());
fAngularDistributionAntiprotonCMSCosTheta->Fill(lorentzvectorantiproton.CosTheta());
fAngularDistributionPionplusCMSCosTheta->Fill(lorentzvectorpionplus.CosTheta());
fAngularDistributionPionminusCMSCosTheta->Fill(lorentzvectorpionminus.CosTheta());

fAngularDistributionLambdaCMSTheta->Fill(lorentzvectorlambda.Theta()*TMath::RadToDeg());
fAngularDistributionLambdabarCMSTheta->Fill(lorentzvectorlambdabar.Theta()*TMath::RadToDeg());
fAngularDistributionProtonCMSTheta->Fill(lorentzvectorproton.Theta()*TMath::RadToDeg());
fAngularDistributionAntiprotonCMSTheta->Fill(lorentzvectorantiproton.Theta()*TMath::RadToDeg());
fAngularDistributionPionplusCMSTheta->Fill(lorentzvectorpionplus.Theta()*TMath::RadToDeg());
fAngularDistributionPionminusCMSTheta->Fill(lorentzvectorpionminus.Theta()*TMath::RadToDeg());

TVector3 momentumvectorlambda = momentumproton+momentumpionminus;
TVector3 momentumvectorlambdabar= momentumantiproton+momentumpionplus;

double protontransversalcomponent = momentumproton.Perp(momentumvectorlambda);
double pionminustransversalcomponent = momentumpionminus.Perp(momentumvectorlambda);

double antiprotontransversalcomponent = momentumantiproton.Perp(momentumvectorlambdabar);
double pionplustransversalcomponent = momentumpionplus.Perp(momentumvectorlambdabar);

double protonlongitudinalcomponent = momentumproton*momentumvectorlambda;
double pionminuslongitudinalcomponent = momentumpionminus*momentumvectorlambda;

double antiprotonlongitudinalcomponent = momentumantiproton*momentumvectorlambdabar;
double pionpluslongitudinalcomponent = momentumpionplus*momentumvectorlambdabar;

alphalambda=((protonlongitudinalcomponent-pionminuslongitudinalcomponent)/(protonlongitudinalcomponent+pionminuslongitudinalcomponent));
alphalambdabar=((pionpluslongitudinalcomponent-antiprotonlongitudinalcomponent)/(antiprotonlongitudinalcomponent+pionpluslongitudinalcomponent));

fArmenterosPodolanskiPlot->Fill(alphalambda,protontransversalcomponent);
fArmenterosPodolanskiPlot->Fill(alphalambdabar,antiprotontransversalcomponent);

momentumlambda = (momentumproton+momentumpionminus).Mag();
momentumlambdabar = (momentumantiproton+momentumpionplus).Mag();

fPositionDecayVertexLambdaLambdaBar2D->Fill(positionprotonz,TMath::Sqrt(positionprotonx*positionprotonx+positionprotony*positionprotony));
fPositionDecayVertexLambdaLambdaBar2D->Fill(positionantiprotonz,TMath::Sqrt(positionantiprotonx*positionantiprotonx+positionantiprotony*positionantiprotony));
fPositionDecayVertexLambda2D->Fill(positionprotonz,TMath::Sqrt(positionprotonx*positionprotonx+positionprotony*positionprotony));
fPositionDecayVertexLambdaBar2D->Fill(positionantiprotonz,TMath::Sqrt(positionantiprotonx*positionantiprotonx+positionantiprotony*positionantiprotony));

if(positionprotonz>positionantiprotonz)
{
	fPositionBothLambdasDecayed->Fill(positionprotonz,TMath::Sqrt(positionprotonx*positionprotonx+positionprotony*positionprotony));
}
else
{
	fPositionBothLambdasDecayed->Fill(positionantiprotonz,TMath::Sqrt(positionantiprotonx*positionantiprotonx+positionantiprotony*positionantiprotony));
}

fMeanLifeLambdaCMS->Fill(positionproton*1.115683/momentumlambda);
fMeanLifeLambdaBarCMS->Fill(positionantiproton*1.115683/momentumlambdabar);
fMeanLifeLambdaLAB->Fill(positionproton);
fMeanLifeLambdaBarLAB->Fill(positionantiproton);

if(fVerbose>0)
{
	std::cout << "============= End PndMCTestLambdaLambdaBar::Exec for Event " << eventcounter << std::endl;
}
eventcounter++;
}

void PndMCTestLambdaLambdaBar::Finish()
{
double NumberOfEventsBeforeDiskPos0 = fPositionBothLambdasDecayed->Integral(fPositionBothLambdasDecayed->GetXaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetXaxis()->FindBin(40.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(20.)) ;
double NumberOfEventsBeforeDiskPos1 = fPositionBothLambdasDecayed->Integral(fPositionBothLambdasDecayed->GetXaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetXaxis()->FindBin(37.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(20.)) ;
double NumberOfEventsBeforeDiskPos2 = fPositionBothLambdasDecayed->Integral(fPositionBothLambdasDecayed->GetXaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetXaxis()->FindBin(57.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(20.)) ;
double NumberOfEventsBeforeDiskPos3 = fPositionBothLambdasDecayed->Integral(fPositionBothLambdasDecayed->GetXaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetXaxis()->FindBin(77.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(0.),fPositionBothLambdasDecayed->GetYaxis()->FindBin(20.)) ;

double PercentOfEventsBeforeDiskPos0 = NumberOfEventsBeforeDiskPos0*100/fPositionBothLambdasDecayed->GetEntries();
double PercentOfEventsBeforeDiskPos1 = NumberOfEventsBeforeDiskPos1*100/fPositionBothLambdasDecayed->GetEntries();
double PercentOfEventsBeforeDiskPos2 = NumberOfEventsBeforeDiskPos2*100/fPositionBothLambdasDecayed->GetEntries();
double PercentOfEventsBeforeDiskPos3 = NumberOfEventsBeforeDiskPos3*100/fPositionBothLambdasDecayed->GetEntries();

double NumberOfHitsBeforeDiskPos0 = fPositionDecayVertexLambdaLambdaBar2D->Integral(fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(40.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(20.)) ;
double NumberOfHitsBeforeDiskPos1 = fPositionDecayVertexLambdaLambdaBar2D->Integral(fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(37.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(20.)) ;
double NumberOfHitsBeforeDiskPos2 = fPositionDecayVertexLambdaLambdaBar2D->Integral(fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(57.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(20.)) ;
double NumberOfHitsBeforeDiskPos3 = fPositionDecayVertexLambdaLambdaBar2D->Integral(fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetXaxis()->FindBin(77.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(0.),fPositionDecayVertexLambdaLambdaBar2D->GetYaxis()->FindBin(20.)) ;

double PercentOfHitsBeforeDiskPos0 = NumberOfHitsBeforeDiskPos0*100/fPositionDecayVertexLambdaLambdaBar2D->GetEntries();
double PercentOfHitsBeforeDiskPos1 = NumberOfHitsBeforeDiskPos1*100/fPositionDecayVertexLambdaLambdaBar2D->GetEntries();
double PercentOfHitsBeforeDiskPos2 = NumberOfHitsBeforeDiskPos2*100/fPositionDecayVertexLambdaLambdaBar2D->GetEntries();
double PercentOfHitsBeforeDiskPos3 = NumberOfHitsBeforeDiskPos3*100/fPositionDecayVertexLambdaLambdaBar2D->GetEntries();


std::cout << "++++++++++++++++++++++++++" << std::endl;
std::cout << "Integral till 40 Hits " <<	PercentOfHitsBeforeDiskPos0 << " %" << std::endl;
std::cout << "Integral till 37 Hits " <<	PercentOfHitsBeforeDiskPos1 << " %" <<std::endl;
std::cout << "Integral till 57 Hits " <<	PercentOfHitsBeforeDiskPos2 << " %" <<std::endl;
std::cout << "Integral till 77 Hits " <<	PercentOfHitsBeforeDiskPos3 << " %" <<std::endl;
std::cout << "++++++++++++++++++++++++++" << std::endl;

std::cout << "++++++++++++++++++++++++++" << std::endl;
std::cout << "Integral till 40 Events " <<	PercentOfEventsBeforeDiskPos0 << " %" << std::endl;
std::cout << "Integral till 37 Events " <<	PercentOfEventsBeforeDiskPos1 << " %" <<std::endl;
std::cout << "Integral till 57 Events " <<	PercentOfEventsBeforeDiskPos2 << " %" <<std::endl;
std::cout << "Integral till 77 Events " <<	PercentOfEventsBeforeDiskPos3 << " %" <<std::endl;
std::cout << "++++++++++++++++++++++++++" << std::endl;


gDirectory->mkdir("MC_Analysis");
gDirectory->cd("MC_Analysis");

	  	gDirectory->mkdir("Angular_Distribution");
	  	gDirectory->cd("Angular_Distribution");

	  	gDirectory->mkdir("CMS");
	  	gDirectory->cd("CMS");

	  	fAngularDistributionLambdaCMSCosTheta->Write();
	  	fAngularDistributionLambdabarCMSCosTheta->Write();
	  	fAngularDistributionProtonCMSCosTheta->Write();
	  	fAngularDistributionAntiprotonCMSCosTheta->Write();
	  	fAngularDistributionPionminusCMSCosTheta->Write();
	  	fAngularDistributionPionplusCMSCosTheta->Write();

	  	fAngularDistributionLambdaCMSTheta->Write();
	  	fAngularDistributionLambdabarCMSTheta->Write();
	  	fAngularDistributionProtonCMSTheta->Write();
	  	fAngularDistributionAntiprotonCMSTheta->Write();
	  	fAngularDistributionPionminusCMSTheta->Write();
	  	fAngularDistributionPionplusCMSTheta->Write();

	  	gDirectory->cd("../");
	  	gDirectory->mkdir("LAB");
	  	gDirectory->cd("LAB");

	  	fAngularDistributionLambdaLABCosTheta->Write();
	  	fAngularDistributionLambdabarLABCosTheta->Write();
	  	fAngularDistributionProtonLABCosTheta->Write();
	  	fAngularDistributionAntiprotonLABCosTheta->Write();
	  	fAngularDistributionPionminusLABCosTheta->Write();
	  	fAngularDistributionPionplusLABCosTheta->Write();

	  	fAngularDistributionLambdaLABTheta->Write();
	  	fAngularDistributionLambdabarLABTheta->Write();
	  	fAngularDistributionProtonLABTheta->Write();
	  	fAngularDistributionAntiprotonLABTheta->Write();
	  	fAngularDistributionPionminusLABTheta->Write();
	  	fAngularDistributionPionplusLABTheta->Write();
	  	gDirectory->cd("../../");


	  	fPzVsPtLambdabar->Write();
	  	fPzVsPtLambda->Write();

	  	fPositionMvdMcPoint->Write();

	    fDetectorCoverageMVDMCPointsProton->Write();
	    fDetectorCoverageMVDMCPointsAntiproton->Write();
	    fDetectorCoverageMVDMCPointsPionplus->Write();
	    fDetectorCoverageMVDMCPointsPionminus->Write();

	  	fMeanLifeLambdaCMS->Write();
	  	fMeanLifeLambdaBarCMS->Write();

	  	fPositionDecayVertexLambdaLambdaBar2D->Write();
	  	fPositionDecayVertexLambdaBar2D->Write();
	  	fPositionDecayVertexLambda2D->Write();
	  	fMeanLifeLambdaLAB->Write();
	  	fMeanLifeLambdaBarLAB->Write();
	  	fPositionBothLambdasDecayed->Write();
	  	fArmenterosPodolanskiPlot->Write();

	  	fMCTracksPerEvent->Write();

	    gDirectory->cd("/");
}

ClassImp(PndMCTestLambdaLambdaBar);