// ************************************************************************
//
// Template for an Analysis Task, 
// 
// for the December 2013 CM Tutorial 
//
// For further info see also
//
// http://panda-wiki.gsi.de/cgi-bin/viewauth/Computing/PandaRootRhoTutorial
// http://panda-wiki.gsi.de/cgi-bin/view/Computing/PandaRootAnalysisJuly13
//
// K.Goetzen 11/2013
//
// ************************************************************************


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

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

// 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"
#include "RhoTuple.h"

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


// -----   Default constructor   -------------------------------------------
PndScrutAnaTask::PndScrutAnaTask(double pbarmom, TString outname) :
  FairTask("Panda Scrutiny Analysis Task") 
{ 
	double mp=0.938272;
	fIni.SetXYZT(1e-12,1e-12,pbarmom, sqrt(pbarmom*pbarmom+mp*mp)+mp);
	hc_rest = fIni.BoostVector();
	// hc_rest.SetX(1e-12); hc_rest.SetY(1e-12);
	std::cout<<"pbarmom = "<<pbarmom<<std::endl;
	fOutName = outname;
	fMethod = "PHSP";
}
// -------------------------------------------------------------------------


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


// -----   Public method Init   --------------------------------------------
InitStatus PndScrutAnaTask::Init() 
{		
	// *** initialize PndAnalysis object
	fAnalysis = new PndAnalysis();
	
	// *** reset the event counter
	fEvtCount = 0;

	// *******
	// ******* PREPARE/CREATE THE STUFF YOU NEED
	// *******
	
	fPdg = TDatabasePDG::Instance();
	
	// ***
	// *** Prepare RhoTuple output  
	// ***
	TDirectory *dir = gDirectory;	
	fFile=new TFile(fOutName,"RECREATE");
	fFile->cd();
	
	// *** create some ntuples
	ntp1 = new RhoTuple("ntp1", "pi0 analysis");
	ntp2 = new RhoTuple("ntp2", "pi+pi- pair analysis");
	ntp5 = new RhoTuple("ntp5", "eta analysis");
	ntp3 = new RhoTuple("ntp3", "h_c analysis");
	ntp4 = new RhoTuple("ntp4", "h_c 4C fit");
	ntp6 = new RhoTuple("ntp6", "pi+ analysis");
	ntp7 = new RhoTuple("ntp7", "pi- analysis");
	//	ntpall = new RhoTuple("ntpall", "all rec info");
	nmc  = new RhoTuple("nmc",  "mctruth info");
	hchi2_4C = new TH1D("hchi2_4C",";#chi^{2}",2e2,0,1e2);

	// assign RhoTuples to outfile
	if (ntp1) ntp1->GetInternalTree()->SetDirectory(gDirectory);
	if (ntp2) ntp2->GetInternalTree()->SetDirectory(gDirectory);
	if (ntp3) ntp3->GetInternalTree()->SetDirectory(gDirectory);
	if (ntp4) ntp4->GetInternalTree()->SetDirectory(gDirectory);
	if (ntp5) ntp5->GetInternalTree()->SetDirectory(gDirectory);
	if (ntp6) ntp6->GetInternalTree()->SetDirectory(gDirectory);
	if (ntp7) ntp7->GetInternalTree()->SetDirectory(gDirectory);
	if (nmc)  nmc->GetInternalTree()->SetDirectory(gDirectory);
	if(hchi2_4C) hchi2_4C->SetDirectory(gDirectory);
	// *** restore original gDirectory
	dir->cd();
	
	return kSUCCESS;
}

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

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


// -----   Public method Exec   --------------------------------------------
void PndScrutAnaTask::Exec(Option_t* opt)
{
  // *** some variables
  int i=0,j=0,k=0;
	
  // *** necessary to read the next event
  fAnalysis->GetEvent();
	
  // *** print event counter
  if (!(++fEvtCount%100)) cout << "evt "<<fEvtCount<<endl;
	
  // *******
  // ******* PUT ANALYSIS CODE HERE
  // *******
	
  TString pidalg = "PidChargedProbability";
  // *** QA tool for simple dumping of analysis results in RhoRuple
  // *** WIKI: https://panda-wiki.gsi.de/foswiki/bin/view/Computing/PandaRootAnalysisJuly13#PndRhoTupleQA
  PndRhoTupleQA qa(fAnalysis,fIni.P());

  // *** RhoCandLists for the analysis
  RhoCandList all, gamma, pi0, pi0pair, eta, piplus, piminus,pichpair, rho, hcfit, hc, allCharged, allNeutral, mclist, hcall, hcreduced, hctmp, pi0tmp, piminustmp, piplustmp, rhotmp, etatmp;
  
  // *** Select with PID info pidalg and ('All'); type and mass are set 		
  fAnalysis->FillList(gamma,  "Neutral");
  fAnalysis->FillList(piplus,  "PionAllPlus",  pidalg);
  fAnalysis->FillList(piminus, "PionAllMinus",  pidalg);
  // fAnalysis->FillList(piplus,  "PionLoosePlus",  pidalg);
  // fAnalysis->FillList(piminus, "PionLooseMinus",  pidalg);
  int ngamma = gamma.GetNumberOfTracks();
  int npipl = piplus.GetNumberOfTracks();
  int npimn = piminus.GetNumberOfTracks();

  // *** Setup event shape object
  fAnalysis->FillList(all,   "All", pidalg);
  PndEventShape evsh(all, fIni, 0.05, 0.1); //PndEventShape(RhoCandList &l, TLorentzVector cms, double neutMinE=0.0, double chrgMinP=0.0);    
  // *** store MC info in ntuple
  fAnalysis->FillList(mclist, "McTruth");
  nmc->Column("ev", (Int_t) fEvtCount);
  qa.qaMcList("",  mclist, nmc);
  nmc->DumpData();
 

   // *** combinatorics for pi0 -> gamma gamma
  pi0.Combine(gamma, gamma);
  pi0.SetType(111);
  
  // *** some rough mass selection on pi0
  // *** Mass selector for the pi0 cands
  double m0_pi0 = fPdg->GetParticle("pi0")->Mass();   // Get nominal PDG mass of the h_c
  //  double pi0MassWind = 0.05;//TEST
  double pi0MassWind = 1e6;//TEST
  RhoMassParticleSelector *pi0MassSel=new RhoMassParticleSelector("pi0",m0_pi0,pi0MassWind);
  pi0.Select(pi0MassSel);
  int npi0 = pi0.GetLength();
  if(ngamma==6 && npipl==1 && npimn==1){
    //    cout<<"npi0 = "<<npi0<<" npipl = "<<npipl<<" npimn = "<<npimn<<endl;
    
    // NEW! ###########################################
    // // *** combinatorics for hc -> (pi+ pi-)3pi0
    pichpair.Combine(piminus,piplus);
    hcall.Combine(piminus,piplus,pi0,pi0,pi0);
    if(hcall.GetLength()>0){
      // *** write ntuple for the hc reconstruction
      for (j=0;j<hcall.GetLength();++j) 
	{
	  if(fMethod=="OmegaEta"){
	    RhoCandidate *pipl = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(0));
	    RhoCandidate *pimn = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(1));
	    RhoCandidate *pi0_1 = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(2));
	    RhoCandidate *pi0_2 = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(3));
	    RhoCandidate *pi0_3 = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(4));
	    pimn->SetType(-211);
	    pipl->SetType(211);
	    pi0_1->SetType(111);
	    pi0_2->SetType(111);
	    pi0_3->SetType(111);
	    RhoCandidate *combCand_eta = pi0_1->Combine(pi0_2,pi0_3);
	    double m2eta = combCand_eta->M()*combCand_eta->M();
	    if(m2eta>0.35 || m2eta<0.25) continue; //only eta
	    eta.Append(combCand_eta);
	    RhoCandidate *combCand_rho = pipl->Combine(pimn);
	    rho.Append(combCand_rho);
	    RhoCandidate *combCand_hc = combCand_rho->Combine(combCand_eta);
	    hcreduced.Append(combCand_hc);
	  }
	  else{
	    hcreduced.Append(hcall[j]);
	  }
	}
      // *** do 4c fit
      std::vector<double> chi2;
      for (j=0;j<hcreduced.GetLength();++j) 
	{
	  //	fIni.Print();
	  Pnd4CFitter fitter(hcreduced[j],fIni);     // *** instantiate the 4C fitter
	  //	fitter.FitConserveMasses();           // *** perform the fit
	  fitter.Fit();           // *** perform the fit with floating of the masses of the electrons and pions
	  double chi2_4c = fitter.GetChi2();    // *** access the chi2 of the fit
	  //	  cout<<"chi2_4c = "<<chi2_4c<<endl;
	  //	fitter.PrintTree(hcreduced[j],10);
	  RhoCandidate *tmpfit = hcreduced[j]->GetFit(); // *** access the fitted candidate            
	  //    if (chi2_4c<20 && tmpfit)               // *** when good enough, fill some histos
	  if (tmpfit){               // *** when good enough, fill some histos
	    hchi2_4C->Fill(chi2_4c);
	    double chi2Wind = 2e1;//TEST
	    //	  double chi2Wind = 1e9;//TEST
	    if(chi2_4c<chi2Wind){
	      hctmp.Add(tmpfit);
	      chi2.push_back(chi2_4c);
	    }
	  }
	}
      // cout<<"hctmp.GetLength()="<<hctmp.GetLength()<<endl;
      if(hctmp.GetLength()>0){
	double ch2_min=20; double j_min=-100;      
	for (j=0;j<hctmp.GetLength();++j){
	  //	cout<<"current ch2 "<<chi2[j]<<endl;
	  if(ch2_min>chi2[j]){
	    ch2_min = chi2[j];
	    j_min = j;
	  }
	}
	// cout<<"cand #"<<j_min<<" has min chi2 = "<<ch2_min<<endl;
	hcfit.Add(hctmp[j_min]);
	TLorentzVector lv_hc = hctmp[j_min]->P4(); 
	hc_rest = lv_hc.BoostVector();
	//  hc_rest.SetX(0.1);  hc_rest.SetY(0.1);
      }
      //      cout<<"hcfit.GetLength()="<<hcfit.GetLength()<<endl;      
      hcfit.SetType(10443);
      int npsimct = fAnalysis->McTruthMatch(hcfit); // match the whole list to count #matches (should be only 1)
      // *** write ntuple for the hc fit
      for (j=0;j<hcfit.GetLength();++j) 
	{
	  // some general info about event (actually written for each candidate!)
	  ntp4->Column("ev",		(Float_t) i);
	  ntp4->Column("cand",	(Float_t) j);
	  ntp4->Column("ncand",   (Float_t) hcfit.GetLength());
	  ntp4->Column("nmct",    (Float_t) npsimct);
	  
	  // store info about initial 4-vector
	  qa.qaP4("beam", fIni, ntp4);
	  
	  // dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
	  qa.qaComp("hcfit", hcfit[j], ntp4);
	  // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
	  RhoCandidate *truth = hcfit[j]->GetMcTruth();
	  TLorentzVector lv;
	  if (truth) lv = truth->P4();
	  qa.qaP4("trhc", lv, ntp4);
	  ntp4->DumpData();
	}
 
      if(hcfit.GetLength()>0){ 
	int nd = hcfit[0]->NDaughters();
	if(nd==5){ //PHSP
	  //store info about final state particles
	  //	  cout<<"pi0_1"<<endl;
	  RhoCandidate *pipl_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(0));
	  piplustmp.Add(pipl_tmp);
	  RhoCandidate *pimn_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(1));
	  piminustmp.Add(pimn_tmp);
	  RhoCandidate *pi0tmp_1 = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(2));
	  pi0tmp.Add(pi0tmp_1);
	  //	cout<<"pi0_2"<<endl;
	  RhoCandidate *pi0tmp_2 = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(3));
	  pi0tmp.Add(pi0tmp_2);
	  //	cout<<"pi0_3"<<endl;
	  RhoCandidate *pi0tmp_3 = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(4));
	  pi0tmp.Add(pi0tmp_3);
	  RhoCandidate *rho_tmp = pipl_tmp->Combine(pimn_tmp);
	  rhotmp.Add(rho_tmp);
	  RhoCandidate *eta_tmp = pi0tmp_1->Combine(pi0tmp_2,pi0tmp_3);
	  etatmp.Add(eta_tmp);
	}
	if(nd==2){ //OmegaEta
	  RhoCandidate *pipl_tmp =RhoFactory::Instance()->NewCandidate((hcfit[0]->Daughter(0))->Daughter(0));
	  RhoCandidate *pimn_tmp = RhoFactory::Instance()->NewCandidate((hcfit[0]->Daughter(0))->Daughter(1));
	  piplustmp.Add(pipl_tmp);
	  piminustmp.Add(pimn_tmp);
	  RhoCandidate *pi0tmp_1 = RhoFactory::Instance()->NewCandidate((hcfit[0]->Daughter(1))->Daughter(0));
	  pi0tmp.Add(pi0tmp_1);
	  RhoCandidate *pi0tmp_2 = RhoFactory::Instance()->NewCandidate((hcfit[0]->Daughter(1))->Daughter(1));
	  pi0tmp.Add(pi0tmp_2);
	  RhoCandidate *pi0tmp_3 = RhoFactory::Instance()->NewCandidate((hcfit[0]->Daughter(1))->Daughter(2));
	  pi0tmp.Add(pi0tmp_3);
	  RhoCandidate *rho_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(0));
	  rhotmp.Add(rho_tmp);
	  RhoCandidate *eta_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(1));
	  etatmp.Add(eta_tmp);
	}
	int npi0mct = fAnalysis->McTruthMatch(pi0tmp); // match the whole list to count #matches (should be only 1)
	// *** write ntuple for the pi0 reconstruction
	for (j=0;j<pi0tmp.GetLength();++j) 
	  {
	    // some general info about event (actually written for each candidate!)
	    ntp1->Column("ev",		(Float_t) i);
	    ntp1->Column("cand",	(Float_t) j);
	    ntp1->Column("ncand",   (Float_t) pi0tmp.GetLength());
	    ntp1->Column("nmct",    (Float_t) npi0mct);
	    
	    // store info about initial 4-vector
	    qa.qaP4("beam", fIni, ntp1);
	    
	    // dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
      //	qa.qaComp("pi0", pi0[j], ntp1);
	    qa.qaPi0("pi0", pi0tmp[j], ntp1);
	    // // dump info about event shapes
	    //		qa.qaEventShapeShort("pi0evsh",&evsh, ntp1);
	    
	    // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
	    RhoCandidate *truth = pi0tmp[j]->GetMcTruth();
	    TLorentzVector lv;
	    if (truth) lv = truth->P4();
	    qa.qaP4("trpi0", lv, ntp1);
	    qa.qaMcDiff("mcdiffpi0", pi0tmp[j], ntp1);
	    TLorentzVector lv_hc = pi0tmp[j]->P4();
	    lv_hc.Boost(-hc_rest);
	    qa.qaP4("pi0_hc_rest_", lv_hc, ntp1);
	    ntp1->DumpData();
	  }
	//fill info about reconstructed pi+/pi-
	int npiplmc = fAnalysis->McTruthMatch(piplustmp); // match the whole list to count #matches (should be only 1)
	//    cout<<"npipl = "<<npipl<<endl;
	for (j=0;j<piplustmp.GetLength();++j) 
	  {
	    ntp6->Column("ev",		(Float_t) i);
	    ntp6->Column("cand",	(Float_t) j);
	    ntp6->Column("ncand",   (Float_t) piplustmp.GetLength());
	    ntp6->Column("nmc",    (Float_t) npiplmc);
	    
	    // store info about initial 4-vector
      qa.qaP4("beam", fIni, ntp6);
      
      // dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
      //	qa.qaComp("pi0", pi0[j], ntp1);
      qa.qaComp("pipl", piplustmp[j], ntp6);
      // // dump info about event shapes
      //		qa.qaEventShapeShort("pi0evsh",&evsh, ntp1);
      
      // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
      RhoCandidate *truth = piplustmp[j]->GetMcTruth();
      TLorentzVector lv;
      if (truth) lv = truth->P4();
      qa.qaP4("trpipl", lv, ntp6);
      qa.qaMcDiff("mcdiffpipl", piplustmp[j], ntp6);
      TLorentzVector lv_hc = piplustmp[j]->P4();
      lv_hc.Boost(-hc_rest);
      qa.qaP4("pipl_hc_rest_", lv_hc, ntp6);
      ntp6->DumpData();
    }
  int npimnmc = fAnalysis->McTruthMatch(piminustmp); // match the whole list to count #matches (should be only 1)
  for (j=0;j<piminustmp.GetLength();++j) 
    {
      ntp7->Column("ev",		(Float_t) i);
      ntp7->Column("cand",	(Float_t) j);
      ntp7->Column("ncand",   (Float_t) piminustmp.GetLength());
      ntp7->Column("nmc",    (Float_t) npimnmc);
      
      // store info about initial 4-vector
      qa.qaP4("beam", fIni, ntp7);
      // dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
      //	qa.qaComp("pi0", pi0[j], ntp1);
      qa.qaComp("pimn", piminustmp[j], ntp7);
      
      // // dump info about event shapes
      //		qa.qaEventShapeShort("pi0evsh",&evsh, ntp1);
      
      // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
      RhoCandidate *truth = piminustmp[j]->GetMcTruth();
      TLorentzVector lv;
      if (truth) lv = truth->P4();
      qa.qaP4("trpimn", lv, ntp7);
      qa.qaMcDiff("mcdiffpimn", piminustmp[j], ntp7);
      TLorentzVector lv_hc = piminustmp[j]->P4();
      lv_hc.Boost(-hc_rest);
      qa.qaP4("pimn_hc_rest_", lv_hc, ntp7);
      ntp7->DumpData();
    }
  //END store info about final state particles


      //Fill ntuples ---------------------------
      // *** write ntuple for the rho reconstruction
      
      //   eta.SetType(221);
      int netat = fAnalysis->McTruthMatch(etatmp); // match the whole list to count #matches (should be only 1)
      //    if(netat!=0) cout<<"netat ="<<netat<<endl;
      for (j=0;j<etatmp.GetLength();++j) 
	{
	  // some general info about event (actually written for each candidate!)
	  ntp5->Column("ev",		(Float_t) i);
	  ntp5->Column("cand",	(Float_t) j);
	  ntp5->Column("ncand",   (Float_t) etatmp.GetLength());
	  ntp5->Column("nmct",    (Float_t) netat);
	  
	  // store info about initial 4-vector
	  qa.qaP4("beam", fIni, ntp5);
	  
	  // dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
	  qa.qaComp("eta", etatmp[j], ntp5);
	  
	  // // dump info about event shapes
	  //	        qa.qaEventShapeShort("etaevsh",&evsh, ntp5);
	  //	cout<<"Daughters of eta ="<<eta[j]->NDaughters()<<endl;
	  // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
	  RhoCandidate *truth = etatmp[j]->GetMcTruth();
	  TLorentzVector lv;
	  if (truth) lv = truth->P4();
	  qa.qaP4("treta", lv, ntp5);
	  // *** store the 4-vector
	  TLorentzVector lv_eta = etatmp[j]->P4();
	  lv_eta.Boost(-hc_rest);
	  qa.qaP4("eta_hc_rest_", lv_eta, ntp5);

	  ntp5->DumpData();
	}
      //   rho.SetType(223); //omega
      rhotmp.SetType(113); //rho
      //   rho.SetType(9010221);//f_0(980)
      int npipairt = fAnalysis->McTruthMatch(rhotmp); // match the whole list to count #matches (should be only 1)
      if(npipairt<1){
	rhotmp.SetType(223); //omega
	npipairt = fAnalysis->McTruthMatch(rhotmp); // match the whole list to count #matches (should be only 1)
	if(npipairt<1){
	  rhotmp.SetType(9010221);//f_0(980)
	  npipairt = fAnalysis->McTruthMatch(rhotmp); // match the whole list to count #matches (should be only 1)
	  if(npipairt<1){
	    rhotmp.SetType(9030221);//f_0(1500)
	    npipairt = fAnalysis->McTruthMatch(rhotmp); // match the whole list to count #matches (should be only 1)
	  }
	}
      }
      //    if(npipairt!=0) cout<<"npipairt ="<<npipairt<<endl;
      for (j=0;j<rhotmp.GetLength();++j) 
	{
	  ntp2->Column("ev",		(Float_t) i);
	  ntp2->Column("cand",	(Float_t) j);
	  ntp2->Column("ncand",   (Float_t) rhotmp.GetLength());
	  ntp2->Column("npipairt",    (Float_t) npipairt);
	  // store info about initial 4-vector
	  qa.qaP4("beam", fIni, ntp2);
	  
	  // dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
	  qa.qaComp("rho", rhotmp[j], ntp2);
	  
	  // // dump info about event shapes
	  //	  qa.qaEventShapeShort("pichpairevsh",&evsh, ntp2);
	  
	  // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
	  RhoCandidate *truth = rhotmp[j]->GetMcTruth();
	  TLorentzVector lv;
	  if (truth) lv = truth->P4();
	  qa.qaP4("trrho", lv, ntp2);
	  
	  // *** store the 4-vector
	  TLorentzVector lv_rho = rhotmp[j]->P4();
	  lv_rho.Boost(-hc_rest);
	  qa.qaP4("rho_hc_rest_", lv_rho, ntp2);
	  
	  ntp2->DumpData();
	}
      
      //*** write ntuple for the eta reconstruction
      //  eta.SetType(221);
      // // *** Mass selector for the eta cands
      // double m0_eta = fPdg->GetParticle(221)->Mass();   // Get nominal PDG mass of the h_c
      // //	double etaMassWind = 1.5;//TEST
      // double etaMassWind = 1e6;//TEST
      // RhoMassParticleSelector *etaMassSel=new RhoMassParticleSelector("eta",m0_eta,etaMassWind);
      // eta.Select(etaMassSel);
      }
      
      // *** write ntuple for the hc reconstruction
      //    hcreduced.SetType(10443);
      hcreduced.SetType(88880);//pbarpsystem
      //    int npsimct = fAnalysis->McTruthMatch(hcreduced); // match the whole list to count #matches (should be only 1)
      for (j=0;j<hcreduced.GetLength();++j){
	// some general info about event (actually written for each candidate!)
	ntp3->Column("ev",		(Float_t) i);
	ntp3->Column("cand",	(Float_t) j);
	ntp3->Column("ncand",   (Float_t) hcreduced.GetLength());
	ntp3->Column("nmct",    (Float_t) npsimct);
	
	// store info about initial 4-vector
	qa.qaP4("beam", fIni, ntp3);
	
	// dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
	qa.qaComp("hc", hcreduced[j], ntp3);
	
	// // dump info about event shapes
	qa.qaEventShapeShort("hcevsh",&evsh, ntp3);
	
	// *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
	RhoCandidate *truth = hcreduced[j]->GetMcTruth();
	TLorentzVector lv;
	if (truth) lv = truth->P4();
	qa.qaP4("trhc", lv, ntp3);
	ntp3->DumpData();
      }
    }
  }
}


// NEW! ###########################################

// // // *** combinatorics for rho0(omega) -> pi+ pi-
// // *** Mass selector for the pi0 cands
// double m0_rho = fPdg->GetParticle(223)->Mass();   // Get nominal PDG mass of the omega
// //	  double rhoMassWind = 0.5;//TEST
// double rhoMassWind = 1e6;//TEST
// RhoMassParticleSelector *rhoMassSel=new RhoMassParticleSelector("omega",m0_rho,rhoMassWind);
// for (j=0;j<piplus.GetLength();++j){
//   for (k=0;k<piminus.GetLength();++k){
//     RhoCandidate *combCand = piplus[j]->Combine(piminus[k]);
//     if(fMethod=="OmegaEta"){
//     double pt = combCand->Pt();
//     double pz = combCand->Pz();
//     double pt_calc = 0.609998 + 0.999846*pz -0.391191*pz*pz + 0.0787386*pz*pz*pz -0.00727785*pz*pz*pz*pz;
//     double diff = fabs(pt-pt_calc);
//     if(diff<0.2)
//       rho.Append(combCand);
//     }
//     else{
//       rho.Append(combCand);
//     }
//   }
// }
// 	// p_perp vs pz
// //width 0.2
// //Chi2                      =      362.215
// //NDf                       =          698
// //p0                        =     0.609998   +/-   0.00210558  
// //p1                        =     0.999846   +/-   0.00636812  
// //p2                        =    -0.391191   +/-   0.00591124  
// //p3                        =    0.0787386   +/-   0.00205302  
// //p4                        =  -0.00727785   +/-   0.000233719 

// //	rho.Combine(piplus, piminus);
// rho.Select(rhoMassSel);
// rho.SetType(223); //omega


// //	pichpair.Combine(piplus, piminus);

// //	pichpair.SetType(113); //rho
// //	pichpair.SetType(223); //omega
//  // *** some rough mass selection on pi0
	  
// //pichpair.Select(rhoMassSel);
//   int npipairt = fAnalysis->McTruthMatch(rho); // match the whole list to count #matches (should be only 1)
//   //cout<<"npipairt = "<<npipairt<<endl;
//   //cout<<"pichpair.GetLength() = "<<rho.GetLength()<<endl;
// // *** write ntuple for the hc reconstruction
// for (j=0;j<rho.GetLength();++j) 
// {
//   ntp2->Column("ev",		(Float_t) i);
//   ntp2->Column("cand",	(Float_t) j);
//   ntp2->Column("ncand",   (Float_t) rho.GetLength());
//   ntp1->Column("npipairt",    (Float_t) npipairt);
//   // store info about initial 4-vector
//   qa.qaP4("beam", fIni, ntp2);
	  
//   // dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
//   qa.qaComp("rho", rho[j], ntp2);
	  
//   // // dump info about event shapes
//   //	  qa.qaEventShapeShort("pichpairevsh",&evsh, ntp2);
	  
//   // *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
//   RhoCandidate *truth = rho[j]->GetMcTruth();
//   TLorentzVector lv;
//   if (truth) lv = truth->P4();
//   qa.qaP4("trrho", lv, ntp2);
//   ntp2->DumpData();
// }

// for (j=0;j<piplus.GetLength();++j){
//   for (k=0;k<piminus.GetLength();++k){
//     for (i=0;i<pi0.GetLength();++i){
//       RhoCandidate *combCand = piplus[j]->Combine(piminus[k],pi0[i]);
//       if(fMethod=="OmegaEta"){
//     double pt = combCand->Pt();
//     double pz = combCand->Pz();
//     double pt_calc = 0.737817 + 0.911204*pz  -0.373268*pz*pz + 0.0790911*pz*pz*pz -0.00771105*pz*pz*pz*pz;
//     double diff = fabs(pt-pt_calc);
//     if(diff<0.2)
//       eta.Append(combCand);
//       }
//       else{
// 	eta.Append(combCand);
//       }
//   }
// }
// }


// // *** combinatorics for eta -> (pi+ pi-) pi0
// // p_perp vs pz
// //width 0.2
// //Chi2                      =      1295.31
// //NDf                       =          685
// //p0                        =     0.737817   +/-   0.000671512 
// //p1                        =     0.911204   +/-   0.00235632  
// //p2                        =    -0.373268   +/-   0.00249091  
// //p3                        =    0.0790911   +/-   0.000942916 
// //p4                        =  -0.00771105   +/-   0.000113898 

// //eta.Combine(piplus, piminus, pi0);

// eta.SetType(221);
// // *** Mass selector for the eta cands
// double m0_eta = fPdg->GetParticle(221)->Mass();   // Get nominal PDG mass of the h_c
// //	double etaMassWind = 1.5;//TEST
// double etaMassWind = 1e6;//TEST
// RhoMassParticleSelector *etaMassSel=new RhoMassParticleSelector("eta",m0_eta,etaMassWind);
// eta.Select(etaMassSel);
// int netat = fAnalysis->McTruthMatch(eta); // match the whole list to count #matches (should be only 1)
// //cout<<"netat ="<<netat<<endl;
// for (j=0;j<eta.GetLength();++j) 
// {
// 	// some general info about event (actually written for each candidate!)
// 	ntp5->Column("ev",		(Float_t) i);
// 	ntp5->Column("cand",	(Float_t) j);
// 	ntp5->Column("ncand",   (Float_t) eta.GetLength());
// 	ntp5->Column("nmct",    (Float_t) netat);
	
// 	// store info about initial 4-vector
// 	qa.qaP4("beam", fIni, ntp5);
			
// 	// dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
// 	qa.qaComp("eta", eta[j], ntp5);
		
// 	// // dump info about event shapes
// 	//	        qa.qaEventShapeShort("etaevsh",&evsh, ntp5);
		
// 	// *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
// 	RhoCandidate *truth = eta[j]->GetMcTruth();
// 	TLorentzVector lv;
// 	if (truth) lv = truth->P4();
// 	qa.qaP4("treta", lv, ntp5);
// 	ntp5->DumpData();
// }

// // *** combinatorics for hc_c -> rho eta
// hc.Combine(rho, eta);		
// //	hc.Combine(pichpairfit, pi0);		
// hc.SetType(10443);

// int npsimct = fAnalysis->McTruthMatch(hc); // match the whole list to count #matches (should be only 1)
// //cout<<"npsimct ="<<npsimct<<endl;
// // *** write ntuple for the hc reconstruction
// for (j=0;j<hc.GetLength();++j) 
// {
// 	// some general info about event (actually written for each candidate!)
// 	ntp3->Column("ev",		(Float_t) i);
// 	ntp3->Column("cand",	(Float_t) j);
// 	ntp3->Column("ncand",   (Float_t) hc.GetLength());
// 	ntp3->Column("nmct",    (Float_t) npsimct);
		
// 	// store info about initial 4-vector
// 	qa.qaP4("beam", fIni, ntp3);
		
// 	// dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
// 	qa.qaComp("hc", hc[j], ntp3);

// 	// // dump info about event shapes
// 	qa.qaEventShapeShort("hcevsh",&evsh, ntp3);
		
// 	// *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
// 	RhoCandidate *truth = hc[j]->GetMcTruth();
// 	TLorentzVector lv;
// 	if (truth) lv = truth->P4();
// 	qa.qaP4("trhc", lv, ntp3);
		
// 	ntp3->DumpData();
// }			

// // *** the lorentz vector of the initial h_c; important for the 4C-fit
// //	TLorentzVector ini(0, 0, 5.6088485, 6.6250557);//TODO: check numbers!
// // ... in event loop ...
// // *** do 4c fit
// for (j=0;j<hc.GetLength();++j) 
//   {
//     Pnd4CFitter fitter(hc[j],fIni);     // *** instantiate the 4C fitter
//     //	    fIni.Print();
//     fitter.FitConserveMasses();           // *** perform the fit
//     //   fitter.Fit();           // *** perform the fit with floating of the masses of the electrons and pions
//     double chi2_4c = fitter.GetChi2();    // *** access the chi2 of the fit
//     //	    cout<<"chi2_4c = "<<chi2_4c<<endl;
//     RhoCandidate *tmpfit = hc[j]->GetFit(); // *** access the fitted candidate            
//     //    if (chi2_4c<20 && tmpfit)               // *** when good enough, fill some histos
//     if (tmpfit){               // *** when good enough, fill some histos
//       hchi2_4C->Fill(chi2_4c);
//       double chi2Wind = 1e1;//TEST
//       //	         double chi2Wind = 1e9;//TEST
//       if(chi2_4c<chi2Wind)
// 	hcfit.Add(tmpfit);
//     }
//   }      
// //	hcfit.SetType(10443);
// //	int nhct = fAnalysis->McTruthMatch(hcfit); // match the whole list to count #matches (should be only 1)
// // cout<<"nhct ="<<nhct<<endl;
// // *** write ntuple for the hc fit
// for (j=0;j<hcfit.GetLength();++j) 
// {
// 	// some general info about event (actually written for each candidate!)
// 	ntp4->Column("ev",		(Float_t) i);
// 	ntp4->Column("cand",	(Float_t) j);
// 	ntp4->Column("ncand",   (Float_t) hcfit.GetLength());
// 	ntp4->Column("nmct",    (Float_t) npsimct);
		
// 	// store info about initial 4-vector
// 	qa.qaP4("beam", fIni, ntp4);
		
// 	// dump information about composite candidate tree recursively (see PndTools/AnalysisTools/PndRhoTupleQA)
// 	qa.qaComp("hcfit", hcfit[j], ntp4);

// 	// // dump info about event shapes
// 	//		qa.qaEventShapeShort("hcfitevsh",&evsh, ntp4);
		
// 	// *** store the 4-vector of the truth matched candidate (or a dummy, if not matched to keep ntuple consistent)
// 	RhoCandidate *truth = hcfit[j]->GetMcTruth();
// 	TLorentzVector lv;
// 	if (truth) lv = truth->P4();
// 	qa.qaP4("trhc", lv, ntp4);
// 	ntp4->DumpData();
// }			
// }

// //Fill ntuple with all info

	
//}


void PndScrutAnaTask::Finish()
{
	
	// *******
	// ******* STORE YOUR HISTOS AND TUPLES
	// *******
	
	fFile->Write();
	fFile->Close();	

}

ClassImp(PndScrutAnaTask)