// ************************************************************************
//
// Code derived from ...
//
// psi(2S) -> J/psi (-> mu+ mu-) pi+ pi- Analysis Example Task
// 
// for the Rho Tutorial, see
// http://panda-wiki.gsi.de/cgi-bin/viewauth/Computing/PandaRootRhoTutorial
//
// K.Goetzen 7/2013
// JGM May 2015
// ************************************************************************


// The header file
#include "PndSolmazAnaTask.h"

// C++ headers
#include <string>
#include <iostream>

// FAIR headers
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"

// ROOT headers
#include "TClonesArray.h"
#include "TVector3.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TDatabasePDG.h"

// RHO headers
#include "RhoCandidate.h"
#include "RhoHistogram/RhoTuple.h"
#include "RhoFactory.h"
#include "RhoMassParticleSelector.h"

// Analysis headers
#include "PndAnalysis.h"
#include "Pnd4CFitter.h"
#include "PndKinVtxFitter.h"
#include "PndKinFitter.h"
#include "PndVtxPoca.h"
#include "PndRhoTupleQA.h"		
#include "PndPidCandidate.h"
		
using std::cout;
using std::endl;


// -----   Default constructor   -------------------------------------------
PndSolmazAnaTask::PndSolmazAnaTask(double pbarmom, TString outname, TString pidalg, TString pidcrit) :
  FairTask("Panda Solmaz Analysis Task") {

  double mp=0.938272;
  fIni.SetXYZT(0,0,pbarmom, sqrt(pbarmom*pbarmom+mp*mp)+mp);

  fOutName = outname;
  fPidAlg = pidalg;
  fPidCrit = pidcrit;
 
}
// -------------------------------------------------------------------------


// -----   Destructor   ----------------------------------------------------
PndSolmazAnaTask::~PndSolmazAnaTask() {
  delete fFile;
 }
// -------------------------------------------------------------------------


// ----- Method to select true PID candidates
int PndSolmazAnaTask::SelectTruePid(PndAnalysis *ana, RhoCandList &l)
{
	int removed = 0;
	
	for (int ii=l.GetLength()-1;ii>=0;--ii)
	{
		if ( !(ana->McTruthMatch(l[ii])) )
		{
			l.Remove(l[ii]);
			removed++;
		}
	}
	
	return removed;
}
// -------------------------------------------------------------------------


// -----   Public method Init   --------------------------------------------
InitStatus PndSolmazAnaTask::Init() 
{		
	// initialize analysis object
	fAnalysis = new PndAnalysis();
	
	// reset the event counter
	fEvtCount = 0;
	
	fPdg = TDatabasePDG::Instance();

	// selectors and all that crap here....

	Double_t m0_pi0 = fPdg->GetParticle("pi0")->Mass();
	pi0MassSel = new RhoMassParticleSelector("pi0",m0_pi0,0.1);

	// output file and ntuple ...

	// *** Prepare RhoTuple output
        // ***
	TDirectory *dir = gDirectory;
        fFile=new TFile(fOutName,"RECREATE");
        fFile->cd();

	hGamMult = new TH1F("hGamMult","Photon multiplicity",100,0,100);
	hPi0Mult = new TH1F("hPi0Mult","Pi0 multiplicity",100,0,100);
	
	nmc = new RhoTuple("nmc",  "mctruth info");
	npi0 = new RhoTuple("npi0", "pi0 analysis");

	if (nmc)  nmc->GetInternalTree()->SetDirectory(gDirectory);
	if (npi0)  npi0->GetInternalTree()->SetDirectory(gDirectory);

	dir->cd();

	return kSUCCESS;
}

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

// -------------------------------------------------------------------------


// -----   Public method Exec   --------------------------------------------
void PndSolmazAnaTask::Exec(Option_t* opt)
{
	// necessary to read the next event
	fAnalysis->GetEvent();
	
	if (!(++fEvtCount%100)) cout << "evt "<<fEvtCount<<endl;
	
	PndRhoTupleQA qa(fAnalysis,fIni.P());

	RhoCandList mclist, gamma, pi0;

	// *** store lists
        fAnalysis->FillList(mclist, "McTruth");
	fAnalysis->FillList(gamma, "Neutral");

        hGamMult->Fill(gamma.GetLength());

	// *** Fill MC information

        nmc->Column("ev", (Int_t) fEvtCount);
        qa.qaMcList("",   mclist, nmc);
        nmc->DumpData();

	// **** do other stuff here ...
	
	// *** Create pi0 candidates from gamma list

	RhoCandidate *g[2];
	pi0.Combine(gamma, gamma);
	pi0.Select(pi0MassSel); // Make preselection on invariant mass 
	pi0.SetType("pi0");

        hPi0Mult->Fill(pi0.GetLength());

	for (int i=0;i<pi0.GetLength();i++)
	  {
	    // *** store basic information of pi0 candidate

	    npi0->Column("ev",    (Float_t) fEvtCount);
	    npi0->Column("cand",  (Float_t) i);
	    npi0->Column("ncand", (Float_t) pi0.GetLength());
	    qa.qaP4("beam", fIni, npi0);

	    // call to "qaPi0" includes MC truth match, EMC information, etc.	    

	    qa.qaPi0("ppi0",(pi0)[i],npi0); 

	    // ....

	    npi0->DumpData();

	  }


	// ....

	// *** Cleanup lists

	mclist.Cleanup();	
	gamma.Cleanup();
	pi0.Cleanup();
}


void PndSolmazAnaTask::Finish()
{	
        hPi0Mult->Write();
        hGamMult->Write();
	nmc->Write();
	npi0->Write();
	fFile->Write();
	fFile->Close();
}

ClassImp(PndSolmazAnaTask)