// ************************************************************************
//
// 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>

// 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-19,1e-19,pbarmom, sqrt(pbarmom*pbarmom+mp*mp)+mp);
	hc_rest = fIni.BoostVector();
	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, piplustmp, piminustmp, pi0tmp, etatmp, rhotmp;
	

	
  // *** 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);

    if(ngamma==2 && npipl==2 && npimn==2){
    //	if(0<1){ //TEST
    
    // NEW! ###########################################
    // // *** combinatorics for hc -> 2(pi+ pi-) pi0 [to avoid double counting between pi+pi- pairs!]
    pichpair.Combine(piminus,piplus);
    hcall.Combine(piminus,piminus,piplus,piplus,pi0);
    hcall.SetType(88880);//pbarpsystem
    for (j=0;j<hcall.GetLength();++j) 
      {
	RhoCandidate *combCand_rho[4]; 	  
	bool flagRho[4];	  bool flagEta[4] ;
	    
	RhoCandidate *pimn0 = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(0));
	RhoCandidate *pimn1 = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(1));
	RhoCandidate *pipl0 = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(2));
	RhoCandidate *pipl1 = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(3));
	RhoCandidate *pi0_cand = RhoFactory::Instance()->NewCandidate(hcall[j]->Daughter(4));
	pimn0->SetType(-211);
	pimn1->SetType(-211);
	pipl0->SetType(211);
	pipl1->SetType(211);
	pi0_cand->SetType(111);
	//	pimn0->PrintOn();
	combCand_rho[0] = pimn0->Combine(pipl0);
	combCand_rho[1] = pimn0->Combine(pipl1);
	combCand_rho[2] = pimn1->Combine(pipl1);
	combCand_rho[3] = pimn1->Combine(pipl0);
	if(fMethod=="OmegaEta"){
	  for(int ic=0;ic<4;ic++){
	    flagRho[ic] = false;
	    flagEta[ic] = false;
	  }
	  ///Approach No.1 Peyrou&mass
      	  for(int ic=0;ic<4;ic++){
      	    double pt = combCand_rho[ic]->Pt();
      	    double pz = combCand_rho[ic]->Pz();
      	    double pt_calc = 0.6007 + 1.017*pz -0.4014*pz*pz + 0.0809*pz*pz*pz -0.007409*pz*pz*pz*pz;
	    //    double w_pt_calc = 0.030258-0.0481854*pz+0.0601155*pz*pz-0.0383109*pz*pz*pz+0.0131429*pz*pz*pz*pz
	    //   -0.00225726*pz*pz*pz*pz*pz+0.000152893*pz*pz*pz*pz*pz*pz; //omega
	    double w_pt_calc = 0.096767-0.107489*pz+0.0901286*pz*pz-0.0438637*pz*pz*pz+0.0123739*pz*pz*pz*pz
	      -0.00184024*pz*pz*pz*pz*pz+0.000112331*pz*pz*pz*pz*pz*pz;//rho
      	    double diff_rho = fabs(pt-pt_calc);
      	    //  if(diff_rho<0.1)
	    //	    if(diff_rho<3*w_pt_calc) //TEST
      	      flagRho[ic] = true;


      	    //   else{
      	      //flagRho[ic] = false;
      	      //**check if it eta
      	      RhoCandidate *combCand_eta = pi0_cand->Combine(combCand_rho[ic]);
      	      pt = combCand_eta->Pt();
      	      pz = combCand_eta->Pz();
      	      pt_calc = 0.7468 + 0.8773*pz  -0.34*pz*pz + 0.06752*pz*pz*pz -0.006423*pz*pz*pz*pz;
      	      w_pt_calc = 0.030258-0.0481854*pz+0.0601155*pz*pz-0.0383109*pz*pz*pz+0.0131429*pz*pz*pz*pz-0.00225726*pz*pz*pz*pz*pz
      		+0.000152893*pz*pz*pz*pz*pz*pz;
      	      double diff_eta = fabs(pt-pt_calc);
	      //	      if(diff_eta<3*w_pt_calc) //TEST
      		flagEta[ic] = true;
      	      //	      cout<<"diff_eta = "<<diff_eta<<" diff_rho = "<<diff_rho<<endl;
      	      // if(diff_eta<0.2 && diff_rho>0.2)
      	      // 	flagEta[ic] = true;
      	      // else{
      	      // 	if(diff_rho<0.2 && diff_eta>0.2)
      	      // 	  flagRho[ic] = true;
      	      // 	else{
      	      // 	  if(diff_rho<0.2 && diff_eta<0.2){
      	      // 	    if(diff_rho<diff_eta)	
      // 	      flagRho[ic] = true;
      	      // 	    else 
      	      // 	      flagEta[ic] = true;
      	      // 	  }
      	      // 	}
      	      // }
      	      //    }

      	  }

      	  for(int ik1=0;ik1<4;ik1++){ //use model assumption for hc candidates
      	    for(int ik2=0;ik2<4;ik2++){ //use model assumption for hc candidates
      	      if(ik1==ik2) continue;
      	      if(flagRho[ik1] && flagEta[ik2]){
      		//	cout<<flagRho[ik1]<<" "<<flagEta[ik2]<<endl;
      	      //	      if(flagRho[ik2] && flagEta[ik1]){
      		//		cout<<"combine pi+pi-  #"<<ik1<<" and #"<<ik2<<endl;
      		//		RhoCandidate *combCand_hc = pi0_cand->Combine(combCand_rho[ik1],combCand_rho[ik2]);
      		//	combCand_rho[ik2]->SetType(422); //omega
      		RhoCandidate *combCand_eta = pi0_cand->Combine(combCand_rho[ik2]->Daughter(0),combCand_rho[ik2]->Daughter(1));
      		//TEST
		//	if(combCand_eta->M()>0.6 || combCand_eta->M()<0.5) continue; //NEW
     		//	if(combCand_rho[ik1]->M()>1.1 || combCand_rho[ik1]->M()<0.45) continue; //NEW: omega+rho
		double m2eta = combCand_eta->M()*combCand_eta->M();
		if(m2eta>0.35 || m2eta<0.25) continue; //NEW: omega+rho
      		combCand_eta->SetType(221);
      		eta.Append(combCand_eta);
      		//combCand_rho[ik1]->SetType(223); //omega TEST
		//	combCand_rho[ik1]->SetType(113); //rho TEST
		rho.Append(combCand_rho[ik1]);
      		RhoCandidate *combCand_hc = combCand_rho[ik1]->Combine(combCand_eta);
      		hcreduced.Append(combCand_hc);
      		//	continue;
      	      }
      	    }
      	  }
	  /// END Approach No.1 Peyrou&mass

	  // ///Approach No.2 Breaking momenta
	  // for(int ik1=0;ik1<4;ik1++){ //use model assumption for hc candidates
	  //   for(int ik2=0;ik2<4;ik2++){ //use model assumption for hc candidates
	  //     if(ik1==ik2) continue;
	  //     double m_pi_pi_2 = (combCand_rho[ik1]->M())*(combCand_rho[ik1]->M());
	  //     if(m_pi_pi_2>0.55 && m_pi_pi_2<0.68){
	  // 	flagRho[ik1] = true;
	  // 	RhoCandidate *combCand_eta = pi0_cand->Combine(combCand_rho[ik2]->Daughter(0),combCand_rho[ik2]->Daughter(1));
	  // 	double m_pi_pi_pi_2 = (combCand_eta->M())*(combCand_eta->M());
	  // 	if(m_pi_pi_pi_2>0.28 && m_pi_pi_pi_2<0.32)
	  // 	  flagEta[ik2] = true;
	  // 	if(flagRho[ik1] && flagEta[ik2]){
	  // 	  combCand_eta->SetType(221);
	  // 	  eta.Append(combCand_eta);
	  // 	  combCand_rho[ik1]->SetType(223); //omega
	  // 	  rho.Append(combCand_rho[ik1]);
	  // 	  RhoCandidate *combCand_hc = combCand_rho[ik1]->Combine(combCand_eta);
	  // 	  hcreduced.Append(combCand_hc);
	  // 	}
	  //     }
	  //   }
	  // }
	  // /// END Approach No.2 Breaking momenta
	}
	else{
	
	      hcreduced.Append(hcall[j]);
	      //// --- Combinotoric test
	      // for(int ic1=0;ic1<4;ic1++){
	      // 	for(int ic2=0;ic2<4;ic2++){
	      // 	  if(ic1==ic2) continue;
	      // 	  RhoCandidate *combCand_eta = pi0_cand->Combine(combCand_rho[ic1]->Daughter(0),combCand_rho[ic1]->Daughter(1));
	      // 	  combCand_eta->SetType(221);
	      // 	  eta.Append(combCand_eta);
	      // 	  rho.Append(combCand_rho[ic2]);
	      // 	  combCand_rho[ic2]->SetType(223); //omega
	      // 	  RhoCandidate *combCand_hc = combCand_rho[ic2]->Combine(combCand_eta);
	      // 	  hcreduced.Append(combCand_hc);
	      // 	}
	      //	  }
	      //// END --- Combinotoric test
	}
      }
    // *** 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)
      //	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();
	}

      if(hcfit.GetLength()>0){
	//Fill ntuples 
	int nd = hcfit[0]->NDaughters();
	if(nd==5){ //PHSP
	  //hcall.Combine(piminus,piminus,piplus,piplus,pi0);
	RhoCandidate *pimn0_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(0));
	RhoCandidate *pimn1_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(1));
	RhoCandidate *pipl0_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(2));
	RhoCandidate *pipl1_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(3));
	RhoCandidate *pi0_cand_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(4));
	pimn0_tmp->SetType(-211);
	pimn1_tmp->SetType(-211);
	pipl0_tmp->SetType(211);
	pipl1_tmp->SetType(211);
	pi0_cand_tmp->SetType(111);
	piplustmp.Add(pipl0_tmp); 	piplustmp.Add(pipl1_tmp);
	piminustmp.Add(pimn0_tmp); 	piminustmp.Add(pimn1_tmp); 
	pi0tmp.Add(pi0_cand_tmp);
	rhotmp.Combine(piminustmp,piplustmp);
	etatmp.Combine(pi0tmp,piminustmp,piplustmp);
	}
	if(nd==2){ //OmegaEta
	  //	  cout<<"Hey!"<<endl;
	  // combCand_hc = combCand_rho[ik1]->Combine(combCand_eta);
	  // RhoCandidate *combCand_eta = pi0_cand->Combine(combCand_rho[ik2]->Daughter(0),combCand_rho[ik2]->Daughter(1));
	  //combCand_rho[0] = pimn0->Combine(pipl0);
	  RhoCandidate *rhoc_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(0));
	  rhotmp.Add(rhoc_tmp);
	  RhoCandidate *etac_tmp = RhoFactory::Instance()->NewCandidate(hcfit[0]->Daughter(1));
	  etatmp.Add(etac_tmp);
	  // int ndeta = etac_tmp->NDaughters();
	  // cout<<"ndeta = "<<ndeta<<endl;
	  // int ndrho = rhoc_tmp->NDaughters();
	  // cout<<"ndrho = "<<ndrho<<endl;

	  RhoCandidate *pimn0_tmp = RhoFactory::Instance()->NewCandidate(rhoc_tmp->Daughter(0));
	  RhoCandidate *pipl0_tmp = RhoFactory::Instance()->NewCandidate(rhoc_tmp->Daughter(1));
	  RhoCandidate *pimn1_tmp = RhoFactory::Instance()->NewCandidate(etac_tmp->Daughter(1));
	  RhoCandidate *pipl1_tmp = RhoFactory::Instance()->NewCandidate(etac_tmp->Daughter(2));
	  RhoCandidate *pi0_cand_tmp = RhoFactory::Instance()->NewCandidate(etac_tmp->Daughter(0));
	  pimn0_tmp->SetType(-211);
	  pimn1_tmp->SetType(-211);
	  pipl0_tmp->SetType(211);
	  pipl1_tmp->SetType(211);
	  pi0_cand_tmp->SetType(111);
	  piplustmp.Add(pipl0_tmp); 	piplustmp.Add(pipl1_tmp);
	  piminustmp.Add(pimn0_tmp); 	piminustmp.Add(pimn1_tmp); 
	  pi0tmp.Add(pi0_cand_tmp);
	  
	}

	//	cout<<"Hey now!"<<endl;
	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("pimn", piminustmp[j], ntp7);
  
	    
	    // *** 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();
	  }

    // *** write ntuple for the rho reconstruction
    int netat = fAnalysis->McTruthMatch(etatmp); // match the whole list to count #matches (should be only 1)
    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);
		
	// *** 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);
	TLorentzVector lv_hc = etatmp[j]->P4();
	lv_hc.Boost(-hc_rest);
	qa.qaP4("eta_hc_rest_", lv_hc, 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)
	}
      }
    }
    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);
	TLorentzVector lv_hc = rhotmp[j]->P4();
	lv_hc.Boost(-hc_rest);
	qa.qaP4("rho_hc_rest_", lv_hc, ntp2);
	ntp2->DumpData();
      }

      }


  

    // *** write ntuple for the hc reconstruction
    //    hcreduced.SetType(10443);
    hcreduced.SetType(88880);//pbarpsystem
    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();
    }


    // // *** do 4c fit
    // 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)
    // 	    hcfit.Add(tmpfit);
    // 	}
    //   }      
    
    // // *** 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();
    //   }
		
  }
}


// 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)