// -------------------------------------------------------------------------
// -----                CbmAnaFlow.h source file             -----
// -----                Created 15/05/12  by Selim         -----
// -------------------------------------------------------------------------
#include <iostream>

#include "CbmAnaFlow.h"

#include "TMath.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TAxis.h"
#include "TRandom.h"
#include "TVector3.h"

#include "FairRootManager.h"
#include "FairRunAna.h"

#include "CbmMCEventHeader.h"
#include "CbmMCTrack.h"
#include "CbmVertex.h"
#include "CbmTrackMatch.h"
#include "CbmStsTrack.h"

#include "CbmKFParticlesFinder.h"
#include "CbmKF.h"
/*#include "CbmKFParticle.h"*/
#include "CbmKFVertex.h"
#include "KFParticle.h"
#include "KFParticleSIMD.h"
#include "KFParticleMatch.h"
#include "CbmL1PFFitter.h"
#include "L1Field.h"
#include "KFMCParticle.h"

#include <iostream>
#include <fstream>
#include <algorithm>
#include <map>
#include <cmath>
#include <vector>

#include "FairTrackParam.h"
#include "CbmKFTrack.h"
#include "CbmKFVertex.h"

#include "TCanvas.h"
#include "TGraphErrors.h"

#include "URun.h"
#include "UEvent.h"
#include "UParticle.h"

using namespace std;
using std::vector;
using std::cout;
using std::endl;

ClassImp(CbmAnaFlow)


// -----   Default constructor   -------------------------------------------
CbmAnaFlow::CbmAnaFlow() : FairTask("EventPlane",1), fHeader(NULL), flistRECOtrack(NULL), flisttrackMATCH(NULL), flistPV(NULL), fRecParticles(NULL), fMCParticles(NULL), fMatchParticles(NULL), flistMCtrack(NULL), fTree_gen(NULL), fEvent_gen(NULL)
{
    fPi = TMath::Pi();
    fmass_proton = 0.938271998;
  
    hPartParam1 = new TH1F ("phi","", 360, -3.2, 3.2);
    hPartParam2 = new TH1F ("mass","", 1000, 0., 10.);
    hPartParam3 = new TH1F ("pt","", 500, 0., 5.);
    hPartParam4 = new TH1F ("Y","", 1000, -5., 5.);
    
    fmode = 1;
    fileName_tree = "";
    fileName_gen = "";
    treeName = "cbmsim_mc";
    fBetaCM = 0.;    
    fGammaCM = 1.;
    fevt_inc = 0;
    
    fdependence = "pt";
    fsel_pdg = 0;
    fsel_factor = "1";
    fsel_Bmin = 0.;
    fsel_Bmax = 0.;
    fy_proj_cm = 0.;
    fy_cm = 0.;
    fEn = 0.;
}

CbmAnaFlow::CbmAnaFlow(const char* name, Int_t verbose, Double_t En) : FairTask(name, verbose), fHeader(NULL), flistRECOtrack(NULL), flisttrackMATCH(NULL), flistPV(NULL), fRecParticles(NULL), fMCParticles(NULL), fMatchParticles(NULL), flistMCtrack(NULL), fTree_gen(NULL), fEvent_gen(NULL)
{
    fPi = TMath::Pi();
    fmass_proton = 0.938271998;
  
    hPartParam1 = new TH1F ("phi","", 360, -3.2, 3.2);
    hPartParam2 = new TH1F ("mass","", 1000, 0., 10.);
    hPartParam3 = new TH1F ("pt","", 500, 0., 5.);
    hPartParam4 = new TH1F ("Y","", 1000, -5., 5.);
    
    fmode = 1;
    fileName_tree = "";
    fileName_gen = "";
    treeName = "cbmsim_mc";
    fBetaCM = 0.;    
    fGammaCM = 1.;
    fevt_inc = 0;
    
    fdependence = "pt";
    fsel_pdg = 0;
    fsel_factor = "1";
    fsel_Bmin = 0.;
    fsel_Bmax = 0.;
    fy_proj_cm = 0.;
    fy_cm = 0.;
    fEn = En;
}

CbmAnaFlow::~CbmAnaFlow()
{
  
  
  //if (fmode == 1)
  //{
  
    if ( fHeader ) 
    {
	fHeader->Delete();
	delete fHeader;
    }    
  
    if ( flistRECOtrack ) 
    {
	flistRECOtrack->Delete();
	delete flistRECOtrack;
    }
    
    if ( flisttrackMATCH ) 
    {
	flisttrackMATCH->Delete();
	delete flisttrackMATCH;
    }
    
    if ( flistPV ) 
    {
	flistPV->Delete();
	delete flistPV;
    }
  
    if ( fRecParticles ) 
    {
	fRecParticles->Delete();
	delete fRecParticles;
    }
    
    if ( fMCParticles ) 
    {
	fMCParticles->Delete();
	delete fMCParticles;
    }
    
    if ( fMatchParticles ) 
    {
	fMatchParticles->Delete();
	delete fMatchParticles;
    }
    
    if ( flistMCtrack ) 
    {
	flistMCtrack->Delete();
	delete flistMCtrack;
    }
    
    if ( fTree_gen )
    {
	fTree_gen->Delete();
	delete fTree_gen;
    }
    
    if ( fEvent_gen )
    {
	fEvent_gen->Delete();
	delete fEvent_gen;
    }
  //}
}

// -----   Public method Init   --------------------------------------------
InitStatus CbmAnaFlow::Init()
{ 
  projRapidityCM();
  
  if (fmode == 0)
  {
    // =========== Get input array
    
    if (fileName_gen == "")
    {
	LOG(FATAL) << "-E- CbmAnaFlow::Init: No generation file!" << FairLogger::endl;
	return kERROR;
    }
    
    cout << "-I- CbmAnaFlow::Init: Opening input file " << fileName_gen << endl;
    TFile *fInputFile = new TFile(fileName_gen);
    if (NULL == fInputFile) 
    {
        LOG(FATAL) << "-E- CbmAnaFlow::Init: Cannot open input file!" << FairLogger::endl;
	return kERROR;	
    }

    // Get run description // -> for elab, plab, etc ... // in fact, CM momenta & energy from UrQMD .. ?
    URun *run = (URun*) fInputFile->Get("run");
    if(NULL == run) 
    {
	LOG(FATAL) << "-E- CbmAnaFlow::Init: No run description in input file!" << FairLogger::endl;
	return kERROR;	
    }
	
    Double_t elab = (TMath::Power(run->GetSqrtS()/run->GetAProj(),2)- 2*TMath::Power(fmass_proton,2))/(2*fmass_proton);
    Double_t plab = TMath::Sqrt(elab*elab - TMath::Power(fmass_proton,2));
    fBetaCM = plab / (elab + fmass_proton);
    fGammaCM = 1./TMath::Sqrt(1. - fBetaCM*fBetaCM);
    
    cout << "-I- CbmAnaFlow::Init: Plab = " << plab << " AGeV, CM boost parameters betaCM = " << fBetaCM << ",  gammaCM = " << fGammaCM << endl;    
    delete run;
    //-> for elab, plab, etc
    
    fTree_gen = (TTree*) fInputFile->Get("events");
    if(NULL == fTree_gen) 
    {
        LOG(FATAL) << "-E- CbmAnaFlow::Init: No event tree in generation input file!" << FairLogger::endl;
	return kERROR;	
    }
    fEvent_gen = new UEvent();
    fTree_gen->SetBranchAddress("event", &fEvent_gen);
    
    //cout << "test0" << endl;
    
    // ============ output Tree
    
    outTree_MC = new TTree("cbmsim_mc","mc particles (generation level)");
      
    outTree_MC->Branch("fphi_to_RP", &fphi_to_RP, "fphi_to_RP/D");  // in fact, from UrQMD, phi_RP = 0. fixed;      
    outTree_MC->Branch("fY", &fY, "fY/D");
    outTree_MC->Branch("fpt", &fpt, "fpt/D");  
    outTree_MC->Branch("fpdg", &fpdg, "fpdf/I");    
    outTree_MC->Branch("fB_MC", &fB_MC, "fB_MC/D");
    
    outTree_MC->SetDirectory(0);
  }       

  if (fmode == 1)
  {
      // =========== Get input array
    
      FairRootManager* ioman = FairRootManager::Instance();
      if ( ! ioman ) 
      {
	cout << "-E- CbmAnaFlow::Init: "
	    << "RootManager not instantised!" << endl;
	return kFATAL;
      }
    
      fHeader = (CbmMCEventHeader*) ioman->GetObject("MCEventHeader.");
      
      if ( !fHeader ) 
      {
	LOG(FATAL) << "-E- CbmAnaFlow::Init: No event header!" << FairLogger::endl;
	return kERROR;
      }
    
      flistMCtrack = (TClonesArray*) ioman->GetObject("MCTrack");
      if ( ! flistMCtrack ) 
      {
	  LOG(FATAL) << "-E- eventPlane::Init: No MC Tracks array!" << FairLogger::endl;
	return kERROR;
      }

      flistRECOtrack = (TClonesArray*) ioman->GetObject("StsTrack");
      if ( ! flistRECOtrack ) 
      {
	  LOG(FATAL) << "-E- eventPlane::Init: No reco. tracks array!" << FairLogger::endl;
	return kERROR;
      }

      flisttrackMATCH = (TClonesArray*) ioman->GetObject("StsTrackMatch");
      if ( ! flisttrackMATCH ) 
      {
	  LOG(FATAL) << "-E- eventPlane::Init: No track match array!" << FairLogger::endl;
	return kERROR;
      }

      flistPV = (CbmVertex*) ioman->GetObject("PrimaryVertex");
      if ( ! flistPV ) 
      {
	  LOG(FATAL) << "-E- eventPlane::Init: No reco. primary vertex array!" << FairLogger::endl;
	return kERROR;
      }
      
      // ============ output Tree  
   
      outTree_MC = new TTree("cbmsim_mc","mc particles (after transport)");
      
      outTree_MC->Branch("fphi_to_RP", &fphi_to_RP, "fphi_to_RP/D");  
      outTree_MC->Branch("fphi_to_EP", &fphi_to_EP, "fphi_to_EP/D");  
      outTree_MC->Branch("fY", &fY, "fY/D");
      outTree_MC->Branch("fpt", &fpt, "fpt/D");  
      outTree_MC->Branch("fpdg", &fpdg, "fpdf/I");
      outTree_MC->Branch("fmass", &fmass, "fmass/D");
      outTree_MC->Branch("fB_MC", &fB_MC, "fB_MC/D");
      outTree_MC->Branch("fmotherID", &fmotherID, "fmotherID/I");
      
      outTree_MC->SetDirectory(0);
      
      outTree_RECO = new TTree("cbmsim_reco","reco. particles");
      
      outTree_RECO->Branch("fphi_to_RP", &fphi_to_RP, "fphi_to_RP/D");  
      outTree_RECO->Branch("fphi_to_EP", &fphi_to_EP, "fphi_to_EP/D");  
      outTree_RECO->Branch("fY", &fY, "fY/D");
      outTree_RECO->Branch("fpt", &fpt, "fpt/D");  
      outTree_RECO->Branch("fpdg", &fpdg, "fpdf/I");
      outTree_RECO->Branch("fmass", &fmass, "fmass/D");
      outTree_RECO->Branch("fB_MC", &fB_MC, "fB_MC/D");
      
      outTree_RECO->SetDirectory(0);
  }

  if (fmode == 2)
  { 
      // =========== Get input array
          
      FairRootManager* ioman = FairRootManager::Instance();
      if ( ! ioman ) 
      {
	cout << "-E- CbmAnaFlow::Init: "
	    << "RootManager not instantised!" << endl;
	return kFATAL;
      }
    
      fHeader = (CbmMCEventHeader*) ioman->GetObject("MCEventHeader.");
      
      if ( !fHeader ) 
      {
	LOG(FATAL) << "-E- CbmAnaFlow::Init: No event header!" << FairLogger::endl;
	return kERROR;
      } 
    
      fRecParticles = dynamic_cast<TClonesArray*>( ioman->GetObject("RecoParticles") );  
      fMCParticles = dynamic_cast<TClonesArray*>( ioman->GetObject("KFMCParticles") );  
      fMatchParticles = dynamic_cast<TClonesArray*>( ioman->GetObject("KFParticleMatch") );  
      flistMCtrack = dynamic_cast<TClonesArray*>( ioman->GetObject("MCTrack") );  

      if ( !fRecParticles || !fMCParticles || !fMatchParticles || !flistMCtrack ) 
      {
	LOG(FATAL) << "-E- CbmAnaFlow::Init: problem with arrays!" << FairLogger::endl;
	return kERROR;
      }     
      
      // ============ output Tree  
   
      outTree_MC = new TTree("cbmsim_mc","mc particles (after transport)");
      
      outTree_MC->Branch("fphi_to_RP", &fphi_to_RP, "fphi_to_RP/D");  
      outTree_MC->Branch("fphi_to_EP", &fphi_to_EP, "fphi_to_EP/D");  
      outTree_MC->Branch("fY", &fY, "fY/D");
      outTree_MC->Branch("fpt", &fpt, "fpt/D");  
      outTree_MC->Branch("fpdg", &fpdg, "fpdf/I");
      outTree_MC->Branch("fmass", &fmass, "fmass/D");
      outTree_MC->Branch("fB_MC", &fB_MC, "fB_MC/D");
      outTree_MC->Branch("fmotherID", &fmotherID, "fmotherID/I");
      
      outTree_MC->SetDirectory(0);
      
      outTree_RECO = new TTree("cbmsim_reco","reco. particles");
      
      outTree_RECO->Branch("fphi_to_RP", &fphi_to_RP, "fphi_to_RP/D");  
      outTree_RECO->Branch("fphi_to_EP", &fphi_to_EP, "fphi_to_EP/D");  
      outTree_RECO->Branch("fY", &fY, "fY/D");
      outTree_RECO->Branch("fpt", &fpt, "fpt/D");  
      outTree_RECO->Branch("fpdg", &fpdg, "fpdf/I");
      outTree_RECO->Branch("fmass", &fmass, "fmass/D");
      outTree_RECO->Branch("fB_MC", &fB_MC, "fB_MC/D");
      
      outTree_RECO->SetDirectory(0);
  }
  
  if (fmode == 3 && fileName_tree == "")
  {
      LOG(FATAL) << "-E- CbmAnaFlow::Init: no input file for histogramming step!" << FairLogger::endl;
      return kERROR;
  }
  
  return kSUCCESS;

}

// -----   Public method Exec   --------------------------------------------
void CbmAnaFlow::Exec(Option_t* opt)
{  
    //if (fmode == 1)
    //{
	//if ( !fHeader ) Fatal("Exec", "No event header");
	//if ( !fRecParticles  || !fMCParticles || !fMatchParticles || !flistMCtrack ) Fatal("Exec", "No particle arrays");
    //}
    
    if (fmode == 0) CreateTreeGen();    
    if (fmode == 1) CreateTreeReco();
    if (fmode == 2) CreateTreeKFPart();
    if (fmode == 3) CreateHisto();
}

void CbmAnaFlow::projRapidityCM()
{
    // Calcul rapidity of c.m. & projectile (in c.m.) -> to select foward Y particles for 1st order event plane
    
    if (fEn == 0.)
    {
      LOG(FATAL) << "-E- CbmAnaFlow::projRapidityCM: Please specify beam kinetic energy" << FairLogger::endl;
      //return kERROR;
    }
    
    const Double_t kProtonMass = 0.938271998;
    Double_t energy_proj = Double_t(fEn) + kProtonMass;
    Double_t mom_proj = sqrt(energy_proj*energy_proj - kProtonMass*kProtonMass);
    Double_t y_proj = 0.5*log((energy_proj+mom_proj)/(energy_proj-mom_proj));
    fy_cm = 0.5*y_proj;
    fy_proj_cm = y_proj - fy_cm;
}

bool CbmAnaFlow::CreateTreeGen()
{ 
   cout << "Number of events in generation file: " << fTree_gen->GetEntries() << endl;  
   
   fTree_gen->GetEntry(fevt_inc);
  
   if(fevt_inc > fTree_gen->GetEntries()) 
   {
     cout << "-I- CbmAnaFlow::CreateTreeGen: End of input file reached" << endl;    
     return kFALSE;
   }
  
   cout << "-I- CbmAnaFlow::CreateTreeGen: Event " << fEvent_gen->GetEventNr() << ", multiplicity " << fEvent_gen->GetNpa() << endl;

   fB_MC = fEvent_gen->GetB();
   cout << "MC B = " << fB_MC << endl; 
   Double_t RP = fEvent_gen->GetPhi();
   cout << "MC RP = " << RP << ", (must be in radian and in [-pi, pi])" << endl;    
   //SELIM: should be =0 in the future (-> random RP in FairPrimaryGenerator)
   // TO CHANGE IF NEEDED
   // RP = RP * fPi/180.;
   // if (RP < -fPi) RP += 2*fPi;
   // if (RP > fPi) RP -= 2*fPi;
   // cout << "MC RP = " << RP << ", (must be in radian and in [-pi, pi])" << endl; 
   
   UParticle *particle;
   Double_t px, py, pz, pz_tmp, energy, energy_tmp;
   
   for (Int_t itrack = 0; itrack < fEvent_gen->GetNpa(); itrack++) 
   {
	// Get particle
	particle = fEvent_gen->GetParticle(itrack);
	
	fpdg = particle->GetPdg();
	
	// momenta & Lorentz boost
	px = particle->Px();
	py = particle->Py();
	fpt = TMath::Sqrt(px*px + py*py); 
	
	fphi_to_RP = TMath::ATan2(py, px);
	fphi_to_RP -= RP;	
	if (fphi_to_RP < -fPi) fphi_to_RP += 2*fPi;
	if (fphi_to_RP > fPi) fphi_to_RP -= 2*fPi;
	
	pz_tmp = particle->Pz();	
	energy_tmp = particle->E();
	
	pz = fGammaCM*(pz_tmp + fBetaCM*energy_tmp);	  // in fact, CM momenta & energy from UrQMD .. ?
	energy = fGammaCM*(energy_tmp + fBetaCM*pz_tmp);
		
	fY = 0.5*TMath::Log( ( energy + pz ) / ( energy - pz ) );
	
	fY -= fy_cm;        // LAB -> CM
	fY /= fy_proj_cm;
	
	// ADD theta -> to set ~same acceptance at generation level & compare with primary reco./mc trak
	
	//cout << "test1 " << endl;
	outTree_MC->Fill();
	//cout << "test2 " << endl;
    }   
    fevt_inc++;
    
    return kTRUE;    
}

bool CbmAnaFlow::CreateTreeReco()
{
    Int_t evtID = fHeader->GetEventID();
    Double_t B = fHeader->GetB();
    Double_t RP = fHeader->GetPhi(); // RP = Reaction plane: MC info
    if ( RP > fPi ) RP -= 2 * fPi;   // en radian!! from SHIELD: [0, 2pi]

    fB_MC = B;
      
    cout << "----- event number: " << evtID << " ------" << endl;
    cout << "----- MC impact parameter: " << B << " ------" << endl;
    cout << "----- MC RP angle: " << RP << ", (must be in radian and in [-pi, pi])" << " ------" << endl;
  
    Int_t trackID, type, motherID;
    Double_t p, px, py, pz, pt, phi, mass, energy, y;
//     Double_t IPx, IPy, IP;
//     Double_t chi2;
//     Int_t NDF;
    CbmStsTrack* track;
    CbmTrackMatch* match;
    CbmMCTrack* mctrack;
    
    FairTrackParam *trackParam;
    TVector3 momRec;
    TVector3 posRec;
    
    // MC tracks
    Int_t nMCtracks = flistMCtrack->GetEntriesFast();
    cout << "evenPlane::CreateTreeReco: # STS MC tracks = " << nMCtracks << endl;
    
    // RECO tracks
    Int_t nSTStracks = flistRECOtrack->GetEntriesFast();
    cout << "evenPlane::CreateTreeReco: # STS reco tracks = " << nSTStracks << endl; 
    
    // Extrapolation track parameters back to primary vertex
    vector<CbmStsTrack> vRTracks;
    vRTracks.resize(nSTStracks);    
    for(int iTr=0; iTr<nSTStracks; iTr++) vRTracks[iTr] = *( (CbmStsTrack*) flistRECOtrack->At(iTr));    
        
    vector<L1FieldRegion> vField;
    vector<float> ChiToPrimVtx;
    CbmL1PFFitter fitter;
    CbmKFVertex kfPV(*flistPV);
    fitter.GetChiToVertex(vRTracks, vField, ChiToPrimVtx, kfPV, 1000000000);    
    
    for (Int_t itrack=0; itrack<nSTStracks; itrack++)
    {
        // =========== RECO. tracks
	
        track = &vRTracks[itrack];
	
	match = (CbmTrackMatch*) flisttrackMATCH->At(itrack);
	trackID = match->GetMCTrackId();

	if (trackID < 0) continue; // ghost track?
	
	//need particle type/mass for rapidity calculation & rotation of pions (anti-v1 flow w.r.t. protons)
	//for now: use MC info for particle ID, later: real particle ID
	//
	mctrack = (CbmMCTrack*) flistMCtrack->At(trackID); 
	type = mctrack->GetPdgCode();
	mass = mctrack->GetMass();      
	
// 	chi2 = track->GetChi2();
// 	NDF = track->GetNDF();
// 	chi2 /= NDF;
	//std::cout << "chi2/NDF = " << chi2 << endl;	
	
	trackParam = track->GetParamFirst();	
	trackParam->Momentum(momRec);
	trackParam->Position(posRec);
	//std::cout << " track Z-StartVertex = " << track->GetParamFirst()->GetZ() << endl;
	
// 	IPx = posRec.X();
// 	IPy = posRec.Y();
// 	IP = TMath::Sqrt(IPx*IPx + IPy*IPy);
	//std::cout << "IPx, IPy, IP = " << IPx << ", " << IPy << ", " << IP << endl;	
	
	px = momRec.X();
	py = momRec.Y();
	pz = momRec.Z();	
	//std::cout << "px, py, pz = " << px << ", " << py << ", " << pz << endl;	
	
	pt = TMath::Sqrt(px*px + py*py);
	p = TMath::Sqrt(px*px + py*py + pz*pz);	
		
	energy = TMath::Sqrt(mass*0.93827203*mass*0.93827203 + p*p);
	y = 0.5 * TMath::Log( ( energy + pz ) / ( energy - pz ) );
	
	phi = TMath::ATan2(py, px);
	//cout << "reco phi = " << phi << ", (must be in radian and in [-pi, pi])" << endl;
	// TO CHANGE IF NEEDED
	//phi = phi * fPi/180.;
	//if (phi < -fPi) phi += 2*fPi;
	//if (phi > fPi) phi -= 2*fPi;
	//cout << "reco phi = " << phi << ", (must be in radian and in [-pi, pi])" << endl;
	
	fphi_to_RP = phi - RP;
	fphi_to_EP = 0.;
	
	if (fphi_to_RP < -fPi) fphi_to_RP += 2*fPi;
	if (fphi_to_RP > fPi) fphi_to_RP -= 2*fPi;
	
	if (fphi_to_EP < -fPi) fphi_to_EP += 2*fPi;
	if (fphi_to_EP > fPi) fphi_to_EP -= 2*fPi;

	fY = y;
	fpt = pt;
	fmass = mass;
	fpdg = type;

	fY -= fy_cm;        // LAB -> CM
	fY /= fy_proj_cm;
	
	outTree_RECO->Fill();
	
	// =========== corresponding MC tracks
	
	y = mctrack->GetRapidity();
	pt = mctrack->GetPt();
	px = mctrack->GetPx();
	py = mctrack->GetPy();
	phi = TMath::ATan2(py, px);
	
	fphi_to_RP = phi - RP;
	fphi_to_EP = 0.;
	
	if (fphi_to_RP < -fPi) fphi_to_RP += 2*fPi;
	if (fphi_to_RP > fPi) fphi_to_RP -= 2*fPi;
	
	if (fphi_to_EP < -fPi) fphi_to_EP += 2*fPi;
	if (fphi_to_EP > fPi) fphi_to_EP -= 2*fPi;
	
	fY = y;
	fpt = pt;
	fmass = mass;
	fpdg = type;
	
	fY -= fy_cm;        // LAB -> CM
	fY /= fy_proj_cm;
	
	outTree_MC->Fill();	
    }
  
  return kTRUE;
}

bool CbmAnaFlow::CreateTreeKFPart()
{
    Int_t evtID = fHeader->GetEventID();
    Double_t B = fHeader->GetB();
    Double_t RP = fHeader->GetPhi(); // RP = Reaction plane: MC info
    if ( RP > fPi ) RP -= 2 * fPi;   // en radian!! from SHIELD: [0, 2pi]

    cout << "----- event number: " << evtID << " ------" << endl;
    cout << "----- MC impact parameter: " << B << " ------" << endl;
    cout << "----- MC RP angle: " << RP << ", (must be in radian and in [-pi, pi])" << " ------" << endl; 
    // SELIM: Should be in [-pi, pi] from FairPrimaryGenerator! now from UnigenGenerator
    // TO CHANGE IF NEEDED
//     RP = RP * fPi/180.;
//     if (RP < -fPi) RP += 2*fPi;
//     if (RP > fPi) RP -= 2*fPi;
//     cout << "----- MC RP angle: " << RP << ", (must be in radian and in [-pi, pi])" << " ------" << endl;
    
    //if (B >= 6. && B < 8.)
    {    
	fB_MC = B;
	
	Double_t M, ErrM;
	//Double_t dL, ErrdL; // decay length
	//Double_t cT, ErrcT; // c*tau
	//Double_t P, ErrP;
	Double_t Pt;
	Double_t Rapidity;
	//Double_t Theta;
	Double_t Phi;
	//Double_t Z;
	
	Double_t px, py;
	
	Int_t nRecoParts = fRecParticles->GetEntriesFast();
	cout << "CbmAnaFlow::CreateTreeKFPart: # reco. part. = " << nRecoParts << endl;

	for (Int_t irec=0; irec<nRecoParts; irec++)
	{
	    // === RECO particle
	    
	    KFParticle partRec = *( (KFParticle*) fRecParticles->At(irec) );
	    KFParticleMatch partMatch = *( (KFParticleMatch*) fMatchParticles->At(irec) );

// 	    if ( ! partRec || ! partMatch ) 
// 	    {
// 	       cout << "CbmAnaFlow::CreateTreeKFPart: no pointer!" << endl;
// 	       continue;
// 	    }
	    
	    if (!partMatch.IsRecoParticle()) continue; // only signals
	    
	    //cout << "Get PDG: " << partRec.GetPDG() << endl; 
	    //partRec.GetMass(M,ErrM);
	    //cout << "Get mass: " << M << endl;
	    
	    int iParticle = fParteff.GetParticleIndex(partRec.GetPDG()); // pdg hypothesis!!
	    if(iParticle < 0) continue;	

	    KFParticle TempPart = partRec;
	    
	    TempPart.GetMass(M,ErrM);
	    //TempPart.GetMomentum(P,ErrP);
	    Pt = TempPart.GetPt();
	    Rapidity = TempPart.GetRapidity();
	    //TempPart.GetDecayLength(dL,ErrdL);
	    //TempPart.GetLifeTime(cT,ErrcT);
	    //Double_t chi2 = TempPart.GetChi2();
	    //Int_t ndf = TempPart.GetNDF();
	    //Double_t prob = TMath::Prob(chi2,ndf);
	    //Theta = TempPart.GetTheta();
	    Phi = TempPart.GetPhi();
	    //Z = TempPart.GetZ();

// 	    KFParticleSIMD tempSIMDPart(TempPart);
// 	    fvec l,dl;
// 	    KFParticleSIMD pv;
// 	    tempSIMDPart.GetDistanceToVertexLine(pv, l, dl);
// 	    dL = sqrt(TempPart.X()*TempPart.X() + TempPart.Y()*TempPart.Y() + TempPart.Z()*TempPart.Z() );	
	    
	    fY = Rapidity;
	    fpt = Pt;
	    fmass = M;
	    fpdg = partRec.GetPDG();
	    
	    cout << "reco Phi = " << Phi << ", (must be in radian and in [-pi, pi])" << endl;
	    // TO CHANGE IF NEEDED
	    //Phi = Phi * fPi/180.;
	    //if (Phi < -fPi) Phi += 2*fPi;
	    //if (Phi > fPi) Phi -= 2*fPi;
	    //cout << "reco Phi = " << Phi << ", (must be in radian and in [-pi, pi])" << endl;
	    
	    fphi_to_RP = Phi - RP;
	    fphi_to_EP = 0.;
	    
	    if (fphi_to_RP < -fPi) fphi_to_RP += 2*fPi;
	    if (fphi_to_RP > fPi) fphi_to_RP -= 2*fPi;
	    
	    if (fphi_to_EP < -fPi) fphi_to_EP += 2*fPi;
	    if (fphi_to_EP > fPi) fphi_to_EP -= 2*fPi;
	    
	    fY -= fy_cm;        // LAB -> CM
	    fY /= fy_proj_cm;
	    
	    outTree_RECO->Fill();
	    
	    hPartParam1->Fill(Phi);
	    hPartParam2->Fill(M);
	    hPartParam3->Fill(Pt);
	    hPartParam4->Fill(Rapidity);
	    
	    // === MC particle	
	    
	    int iMCPart = partMatch.GetMatch();
	    KFMCParticle mcPart = *( (KFMCParticle*) fMCParticles->At(iMCPart) );
	    
	    int iMCTrack = mcPart.GetMCTrackID();
	    CbmMCTrack &mcTrack = *(static_cast<CbmMCTrack*>(flistMCtrack->At(iMCTrack)));
	    int mcDaughterId = mcPart.GetDaughterIds()[0];
	    CbmMCTrack &mcDaughter = *(static_cast<CbmMCTrack*>(flistMCtrack->At(mcDaughterId)));
	    
	    fmotherID = mcTrack.GetMotherId();
	    
	    fY = mcTrack.GetRapidity();
	    fpt = mcTrack.GetPt();
	    fmass = mcTrack.GetMass();
	    fpdg = mcTrack.GetPdgCode();
	    
	    px = mcTrack.GetPx();
	    py = mcTrack.GetPy();
	    Phi = TMath::ATan2(py, px);	    
	    cout << "mc Phi = " << Phi << endl;
	    
	    fphi_to_RP = Phi - RP;		    
	    fphi_to_EP = 0.;
	    
	    if (fphi_to_RP < -fPi) fphi_to_RP += 2*fPi;
	    if (fphi_to_RP > fPi) fphi_to_RP -= 2*fPi;
	    
	    if (fphi_to_EP < -fPi) fphi_to_EP += 2*fPi;
	    if (fphi_to_EP > fPi) fphi_to_EP -= 2*fPi;
	    
	    fY -= fy_cm;        // LAB -> CM
	    fY /= fy_proj_cm;
	    
	    // ADD theta & motherID -> to set ~same acceptance at generation level & compare with primary reco./mc traks  
	    
	    outTree_MC->Fill();
	}
    }  
    
    return kTRUE;
}

void CbmAnaFlow::CreateHisto()
{ 
  cout << "-- file containing Tree for flow analysis: " << fileName_tree << endl;
    
  // later: break run if no file nor tree
  
  TFile *f = new TFile(fileName_tree,"R");
  if (!f) cout << "-E- CbmAnaFlow::CreateHisto: pb with file containing Tree -----" << endl;
  
  TTree* tree = (TTree* ) f->Get(treeName);  
  if (!tree) cout << "-E- CbmAnaFlow::CreateHisto: pb with Tree in file -----" << endl;
  
  TString spdg, ssel_Bmin, ssel_Bmax;
  
  spdg = "";
  ssel_Bmin = "";
  ssel_Bmax = "";
  
  spdg += fsel_pdg;
  ssel_Bmin += fsel_Bmin;
  ssel_Bmax += fsel_Bmax;
  
  if (fdependence == "pt")
  { 
      cout << "========= pt dependence of v" << fsel_factor << " of '" << spdg << "' particles" << endl;
    
      Int_t n = 5;
      Float_t sigma0[5], esigma0[5];
      Float_t smear0[5], esmear0[5];
      
      //TString Yrange = " && fY > -0.5 && fY < 0.5"; // for v2
      TString Yrange = " && (fY < -0.5 || fY > 0.5)"; // for v1
      
      TString cut = "fpdg == " + spdg + " && fB_MC >= " + ssel_Bmin + " && fB_MC < " + ssel_Bmax;
      cut += " && fpt < 0.5" + Yrange;
      cout << "cut: " << cut << endl;
      //TO ADD: cut += motherID ==-1 && theta in 2.5, 25 deg. when using mc from reco track   
      
      tree->Draw("TMath::Cos(" + fsel_factor + "*fphi_to_RP) >> h1", cut);
      TH1F* h1 = (TH1F*) tree->GetHistogram();
      sigma0[0] = h1->GetMean();
      esigma0[0] = h1->GetMeanError();
      h1->Clear();
      
      tree->Draw("fpt >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      smear0[0] = h1->GetMean();
      esmear0[0] = h1->GetMeanError();
      h1->Clear();
      
      cut = "fpdg == " + spdg + " && fB_MC >= " + ssel_Bmin + " && fB_MC < " + ssel_Bmax;
      cut += " && fpt >= 0.5 && fpt < 1." + Yrange;
      cout << "cut: " << cut << endl;
      
      tree->Draw("TMath::Cos(" + fsel_factor + "*fphi_to_RP) >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      sigma0[1] = h1->GetMean();
      esigma0[1] = h1->GetMeanError();
      h1->Clear();
      
      tree->Draw("fpt >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      smear0[1] = h1->GetMean();
      esmear0[1] = h1->GetMeanError();
      h1->Clear();
      
      cut = "fpdg == " + spdg + " && fB_MC >= " + ssel_Bmin + " && fB_MC < " + ssel_Bmax;
      cut += " && fpt >= 1. && fpt < 1.5" + Yrange;
      cout << "cut: " << cut << endl;
      
      tree->Draw("TMath::Cos(" + fsel_factor + "*fphi_to_RP) >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      sigma0[2] = h1->GetMean();
      esigma0[2] = h1->GetMeanError();
      h1->Clear();
      
      tree->Draw("fpt >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      smear0[2] = h1->GetMean();
      esmear0[2] = h1->GetMeanError();
      h1->Clear();
      
      cut = "fpdg == " + spdg + " && fB_MC >= " + ssel_Bmin + " && fB_MC < " + ssel_Bmax;
      cut += " && fpt >= 1.5 && fpt < 2." + Yrange;
      cout << "cut: " << cut << endl;
      
      tree->Draw("TMath::Cos(" + fsel_factor + "*fphi_to_RP) >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      sigma0[3] = h1->GetMean();
      esigma0[3] = h1->GetMeanError();
      h1->Clear();
      
      tree->Draw("fpt >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      smear0[3] = h1->GetMean();
      esmear0[3] = h1->GetMeanError();
      h1->Clear();
      
      cut = "fpdg == " + spdg + " && fB_MC >= " + ssel_Bmin + " && fB_MC < " + ssel_Bmax;
      cut += " && fpt >= 2." + Yrange;
      cout << "cut: " << cut << endl;
      
      tree->Draw("TMath::Cos(" + fsel_factor + "*fphi_to_RP) >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      sigma0[4] = h1->GetMean();
      esigma0[4] = h1->GetMeanError();
      h1->Clear();
      
      tree->Draw("fpt >> h1", cut);
      h1 = (TH1F*) tree->GetHistogram();
      smear0[4] = h1->GetMean();
      esmear0[4] = h1->GetMeanError();
      h1->Clear();
      
      //cout << "test" << endl;

      for (Int_t igraph=0; igraph<5; igraph++) cout << "<pt>, v" + fsel_factor + " = " << smear0[igraph] << ", " << sigma0[igraph] << endl;
      
      TCanvas *c1 = new TCanvas("c1"," sigmaRP vs b ",200,10,700,500);
      c1->SetFillColor(0);
      c1->SetGrid();

      // draw a frame to define the range
      TH1F *hr = c1->DrawFrame(0.,-1.,3.,1.);

      hr->SetXTitle("Tranverse momentum p_{T} [Gev/c]");
      hr->GetXaxis()->CenterTitle();

      hr->SetYTitle("v_{" + fsel_factor + "}");
      hr->GetYaxis()->CenterTitle();

      TGraphErrors *gr0 = new TGraphErrors(n,smear0,sigma0,esmear0,esigma0);
      gr0->SetMarkerColor(kBlack);
      gr0->SetMarkerStyle(25);
      gr0->SetLineColor(kBlack);
      //gr0->Fit("pol1");
      //gr0->Draw("LP");
      gr0->Draw("P");
  }
  
  if (fdependence == "rapidity")
  {
    cout << "========= Rapidity dependence of v" << fsel_factor << " of '" << spdg << "' particles" << endl;
    
    Int_t n = 10;

    Float_t sigma0[10], esigma0[10];
    Float_t smear0[10], esmear0[10];

    TString cut = "&& fpdg == " + spdg + " && fB_MC >= " + ssel_Bmin + " && fB_MC < " + ssel_Bmax;    
    //TO ADD: cut += motherID ==-1 && theta in 2.5, 25 deg. when using mc from reco track   

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h1","fY<-0.8" + cut);
    TH1F* h1 = (TH1F*) tree->GetHistogram();
    sigma0[0] = h1->GetMean();
    esigma0[0] = h1->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h2","fY>=-0.8 && fY<-0.6" + cut);
    TH1F* h2 = (TH1F*) tree->GetHistogram();
    sigma0[1] = h2->GetMean();
    esigma0[1] = h2->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h3","fY>=-0.6 && fY<-0.4" + cut);
    TH1F* h3 = (TH1F*) tree->GetHistogram();
    sigma0[2] = h3->GetMean();
    esigma0[2] = h3->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h4","fY>=-0.4 && fY<-0.2" + cut);
    TH1F* h4 = (TH1F*) tree->GetHistogram();
    sigma0[3] = h4->GetMean();
    esigma0[3] = h4->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h5","fY>=-0.2 && fY<0." + cut);
    TH1F* h5 = (TH1F*) tree->GetHistogram();
    sigma0[4] = h5->GetMean();
    esigma0[4] = h5->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h6","fY>=0. && fY<0.2" + cut);
    TH1F* h6 = (TH1F*) tree->GetHistogram();
    sigma0[5] = h6->GetMean();
    esigma0[5] = h6->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h7","fY>=0.2 && fY<0.4" + cut);
    TH1F* h7 = (TH1F*) tree->GetHistogram();
    sigma0[6] = h7->GetMean();
    esigma0[6] = h7->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h8","fY>=0.4 && fY<0.6" + cut);
    TH1F* h8 = (TH1F*) tree->GetHistogram();
    sigma0[7] = h8->GetMean();
    esigma0[7] = h8->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h9","fY>=0.6 && fY<0.8" + cut);
    TH1F* h9 = (TH1F*) tree->GetHistogram();
    sigma0[8] = h9->GetMean();
    esigma0[8] = h9->GetMeanError();

    tree->Draw("TMath::Cos(" + fsel_factor + " * fphi_to_RP) >>h10","fY>=0.8" + cut);
    TH1F* h10 = (TH1F*) tree->GetHistogram();
    sigma0[9] = h10->GetMean();
    esigma0[9] = h10->GetMeanError();

    //------------
    tree->Draw("fY >>h21","fY<-0.8" + cut);
    TH1F* h21 = (TH1F*) tree->GetHistogram();
    smear0[0] = h21->GetMean();
    esmear0[0] = h21->GetMeanError();

    tree->Draw("fY >>h22","fY>=-0.8 && fY<-0.6" + cut);
    TH1F* h22 = (TH1F*) tree->GetHistogram();
    smear0[1] = h22->GetMean();
    esmear0[1] = h22->GetMeanError();

    tree->Draw("fY >>h23","fY>=-0.6 && fY<-0.4" + cut);
    TH1F* h23 = (TH1F*) tree->GetHistogram();
    smear0[2] = h23->GetMean();
    esmear0[2] = h23->GetMeanError();

    tree->Draw("fY >>h24","fY>=-0.4 && fY<-0.2" + cut);
    TH1F* h24 = (TH1F*) tree->GetHistogram();
    smear0[3] = h24->GetMean();
    esmear0[3] = h24->GetMeanError();

    tree->Draw("fY >>h25","fY>=-0.2 && fY<0." + cut);
    TH1F* h25 = (TH1F*) tree->GetHistogram();
    smear0[4] = h25->GetMean();
    esmear0[4] = h25->GetMeanError();

    tree->Draw("fY >>h26","fY>=0. && fY<0.2" + cut);
    TH1F* h26 = (TH1F*) tree->GetHistogram();
    smear0[5] = h26->GetMean();
    esmear0[5] = h26->GetMeanError();

    tree->Draw("fY >>h27","fY>=0.2 && fY<0.4" + cut);
    TH1F* h27 = (TH1F*) tree->GetHistogram();
    smear0[6] = h27->GetMean();
    esmear0[6] = h27->GetMeanError();

    tree->Draw("fY >>h28","fY>=0.4 && fY<0.6" + cut);
    TH1F* h28 = (TH1F*) tree->GetHistogram();
    smear0[7] = h28->GetMean();
    esmear0[7] = h28->GetMeanError();

    tree->Draw("fY >>h29","fY>=0.6 && fY<0.8" + cut);
    TH1F* h29 = (TH1F*) tree->GetHistogram();
    smear0[8] = h29->GetMean();
    esmear0[8] = h29->GetMeanError();

    tree->Draw("fY >>h210","fY>=0.8" + cut);
    TH1F* h210 = (TH1F*) tree->GetHistogram();
    smear0[9] = h210->GetMean();
    esmear0[9] = h210->GetMeanError();
    
    TCanvas *c1 = new TCanvas("c1"," sigmaRP vs b ",200,10,700,500);
    c1->SetFillColor(0);
    c1->SetGrid();

    for (Int_t igraph=0; igraph<10; igraph++) cout << "<Y>, v" + fsel_factor + " = " << smear0[igraph] << ", " << sigma0[igraph] << endl;
    
    // draw a frame to define the range
    TH1F *hr = c1->DrawFrame(-1.5,-0.3,1.5,0.3);

    hr->SetXTitle("Rapidity (in CM frame, norm. to projectile Y)");
    hr->GetXaxis()->CenterTitle();

    hr->SetYTitle("v_{" + fsel_factor + "}");
    hr->GetYaxis()->CenterTitle();

    TGraphErrors *gr0 = new TGraphErrors(n,smear0,sigma0,esmear0,esigma0);
    gr0->SetMarkerColor(kBlack);
    gr0->SetMarkerStyle(25);
    gr0->SetLineColor(kBlack);
    //gr0->Fit("pol1");
    //gr0->Draw("LP");
    gr0->Draw("P");
  }
  
}

void CbmAnaFlow::Write()
{
  if (fmode == 0)
  {
    /*hPartParam1->Write();    
    hPartParam1->Delete();
    
    hPartParam2->Write();    
    hPartParam2->Delete();
    
    hPartParam3->Write();    
    hPartParam3->Delete();
    
    hPartParam4->Write();    
    hPartParam4->Delete(); */   
    
    outTree_MC->Write("",TObject::kOverwrite);    
  }
  
  if (fmode == 1)
  {
    outTree_MC->Write("",TObject::kOverwrite);
    outTree_RECO->Write("",TObject::kOverwrite);
  }
  
  if (fmode == 2)
  {
    hPartParam1->Write();    
    hPartParam1->Delete();
    
    hPartParam2->Write();    
    hPartParam2->Delete();
    
    hPartParam3->Write();    
    hPartParam3->Delete();
    
    hPartParam4->Write();    
    hPartParam4->Delete();    
    
    outTree_MC->Write("",TObject::kOverwrite);
    outTree_RECO->Write("",TObject::kOverwrite);
  }
}

void CbmAnaFlow::Finish()
{
    // Finish of the task execution
    Write();
}