// ------------------------------------------------------------------
// -----                     CbmHadronAnalysis                  -----
// -----              Created 14/12/2012 by nh                  -----
// ------------------------------------------------------------------
#include <iostream>
using namespace std;

#include "TClonesArray.h"
#include "TH1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TString.h"
#include "TFile.h"
#include "TMath.h"
#include "TRandom.h"
#include "TROOT.h"

#include "FairRootManager.h"
#include "CbmMCEventHeader.h"
#include "CbmMCTrack.h"
#include "CbmStsTrack.h"
#include "CbmStsHit.h"
#include "CbmStsPoint.h"
#include "CbmTofPoint.h"
#include "CbmTofHit.h"
#include "TofData/CbmTofCell.h"
#include "TofTools/CbmTofGeoHandler.h"
#include "CbmTofDetectorId_v12b.h"
#include "CbmGlobalTrack.h"
#include "CbmHadron.h"
#include "CbmHadronAnalysis.h"

#define M2PI 0.019479835
#define M2KA 0.24371698
#define M2PROT 0.88035435
#define clight 29.9792;

//___________________________________________________________________
//
// CbmHadronAnalysis
//
// Task for analysis of hadron spectra
//
// ------------------------------------------------------------------
CbmHadronAnalysis::CbmHadronAnalysis()
  : FairTask("HadronAnalysis"),
    fEvents(0),
    fBeamMomentum(10),
    fPdfFileName(""),
    fMCTracks(NULL),
    fTofPoints(NULL),
    fTofHits(NULL),
    fHadrons(NULL)
{
  CreateHistogramms();
  cout << "CbmHadronAnalysis: Task started " << endl;
}
// ------------------------------------------------------------------

// ------------------------------------------------------------------
CbmHadronAnalysis::CbmHadronAnalysis(const char* name, Int_t verbose)
  : FairTask(name, verbose),
    fEvents(0),
    fBeamMomentum(10),
    fPdfFileName(""),
    fMCTracks(NULL),
    fTofPoints(NULL),
    fTofHits(NULL),
    fHadrons(NULL)
{
    CreateHistogramms();
}
// ------------------------------------------------------------------


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


// ------------------------------------------------------------------
void CbmHadronAnalysis::CreateHistogramms()
{
  // Create histogramms

  Float_t ymin=-1.;
  Float_t ymax=4.;
  Float_t ptmmax=2.5;
  Int_t   ptm_nbx=30;
  Int_t   ptm_nby=30;

  TDirectory * oldir = gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
  gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !

  // generator level

  fa_ptm_rap_gen_pip = new TH2F("ptm_rap_gen_pip","MCTrack-gen pi-plus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_pim = new TH2F("ptm_rap_gen_pim","MCTrack-gen pi-minus;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_kp  = new TH2F("ptm_rap_gen_kp", "MCTrack-gen k-plus;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_km  = new TH2F("ptm_rap_gen_km", "MCTrack-gen k-minus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_p   = new TH2F("ptm_rap_gen_p",  "MCTrack-gen proton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_pbar= new TH2F("ptm_rap_gen_pbar","MCTrack-gen antiproton;y;p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_d   = new TH2F("ptm_rap_gen_d",  "MCTrack-gen deuteron;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_t   = new TH2F("ptm_rap_gen_t",  "MCTrack-gen triton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_h   = new TH2F("ptm_rap_gen_h",  "MCTrack-gen 3he;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_a   = new TH2F("ptm_rap_gen_a",  "MCTrack-gen alpha;   y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_gen_imf = new TH2F("ptm_rap_gen_imf","MCTrack-gen imf;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);

  Float_t v1_nbx=20.;
  Float_t v1_nby=20.;
  Float_t yvmax=1.3;

  fa_v1_rap_gen_pip = new TH2F("v1_rap_gen_pip","MCTrack-gen pi-plus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_pim = new TH2F("v1_rap_gen_pim","MCTrack-gen pi-minus;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_kp  = new TH2F("v1_rap_gen_kp", "MCTrack-gen k-plus;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_km  = new TH2F("v1_rap_gen_km", "MCTrack-gen k-minus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_p   = new TH2F("v1_rap_gen_p",  "MCTrack-gen proton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_pbar= new TH2F("v1_rap_gen_pbar","MCTrack-gen antiproton;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_d   = new TH2F("v1_rap_gen_d",  "MCTrack-gen deuteron;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_t   = new TH2F("v1_rap_gen_t",  "MCTrack-gen triton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_h   = new TH2F("v1_rap_gen_h",  "MCTrack-gen 3he;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_a   = new TH2F("v1_rap_gen_a",  "MCTrack-gen alpha;   y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_gen_imf = new TH2F("v1_rap_gen_imf","MCTrack-gen imf;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  fa_v2_rap_gen_pip = new TH2F("v2_rap_gen_pip","MCTrack-gen pi-plus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_pim = new TH2F("v2_rap_gen_pim","MCTrack-gen pi-minus;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_kp  = new TH2F("v2_rap_gen_kp", "MCTrack-gen k-plus;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_km  = new TH2F("v2_rap_gen_km", "MCTrack-gen k-minus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_p   = new TH2F("v2_rap_gen_p",  "MCTrack-gen proton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_pbar= new TH2F("v2_rap_gen_pbar","MCTrack-gen antiproton;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_d   = new TH2F("v2_rap_gen_d",  "MCTrack-gen deuteron;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_t   = new TH2F("v2_rap_gen_t",  "MCTrack-gen triton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_h   = new TH2F("v2_rap_gen_h",  "MCTrack-gen 3he;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_a   = new TH2F("v2_rap_gen_a",  "MCTrack-gen alpha;   y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_gen_imf = new TH2F("v2_rap_gen_imf","MCTrack-gen imf;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  // TofPoint level

  fa_ptm_rap_poi_pip = new TH2F("ptm_rap_poi_pip","MCTrack-poi pi-plus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_pim = new TH2F("ptm_rap_poi_pim","MCTrack-poi pi-minus;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_kp  = new TH2F("ptm_rap_poi_kp", "MCTrack-poi k-plus;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_km  = new TH2F("ptm_rap_poi_km", "MCTrack-poi k-minus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_p   = new TH2F("ptm_rap_poi_p",  "MCTrack-poi proton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_pbar= new TH2F("ptm_rap_poi_pbar","MCTrack-poi antiproton;y;p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_d   = new TH2F("ptm_rap_poi_d",  "MCTrack-poi deuteron;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_t   = new TH2F("ptm_rap_poi_t",  "MCTrack-poi triton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_h   = new TH2F("ptm_rap_poi_h",  "MCTrack-poi 3he;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_a   = new TH2F("ptm_rap_poi_a",  "MCTrack-poi alpha;   y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_poi_imf = new TH2F("ptm_rap_poi_imf", "MCTrack-poi imf;    y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);

  fa_v1_rap_poi_pip = new TH2F("v1_rap_poi_pip","MCTrack-poi pi-plus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_pim = new TH2F("v1_rap_poi_pim","MCTrack-poi pi-minus;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_kp  = new TH2F("v1_rap_poi_kp", "MCTrack-poi k-plus;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_km  = new TH2F("v1_rap_poi_km", "MCTrack-poi k-minus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_p   = new TH2F("v1_rap_poi_p",  "MCTrack-poi proton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_pbar= new TH2F("v1_rap_poi_pbar","MCTrack-poi antiproton;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_d   = new TH2F("v1_rap_poi_d",  "MCTrack-poi deuteron;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_t   = new TH2F("v1_rap_poi_t",  "MCTrack-poi triton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_h   = new TH2F("v1_rap_poi_h",  "MCTrack-poi 3he;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_a   = new TH2F("v1_rap_poi_a",  "MCTrack-poi alpha;   y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_poi_imf = new TH2F("v1_rap_poi_imf","MCTrack-poi imf;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  fa_v2_rap_poi_pip = new TH2F("v2_rap_poi_pip","MCTrack-poi pi-plus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_pim = new TH2F("v2_rap_poi_pim","MCTrack-poi pi-minus;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_kp  = new TH2F("v2_rap_poi_kp", "MCTrack-poi k-plus;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_km  = new TH2F("v2_rap_poi_km", "MCTrack-poi k-minus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_p   = new TH2F("v2_rap_poi_p",  "MCTrack-poi proton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_pbar= new TH2F("v2_rap_poi_pbar","MCTrack-poi antiproton;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_d   = new TH2F("v2_rap_poi_d",  "MCTrack-poi deuteron;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_t   = new TH2F("v2_rap_poi_t",  "MCTrack-poi triton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_h   = new TH2F("v2_rap_poi_h",  "MCTrack-poi 3he;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_a   = new TH2F("v2_rap_poi_a",  "MCTrack-poi alpha;   y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_poi_imf = new TH2F("v2_rap_poi_imf","MCTrack-poi imf;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  // TofHit level

  fa_ptm_rap_hit_pip = new TH2F("ptm_rap_hit_pip","MCTrack-hit pi-plus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_pim = new TH2F("ptm_rap_hit_pim","MCTrack-hit pi-minus;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_kp  = new TH2F("ptm_rap_hit_kp", "MCTrack-hit k-plus;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_km  = new TH2F("ptm_rap_hit_km", "MCTrack-hit k-minus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_p   = new TH2F("ptm_rap_hit_p",  "MCTrack-hit proton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_pbar= new TH2F("ptm_rap_hit_pbar","MCTrack-hit antiproton;y;p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_d   = new TH2F("ptm_rap_hit_d",  "MCTrack-hit deuteron;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_t   = new TH2F("ptm_rap_hit_t",  "MCTrack-hit triton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_h   = new TH2F("ptm_rap_hit_h",  "MCTrack-hit 3he;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_a   = new TH2F("ptm_rap_hit_a",  "MCTrack-hit alpha;   y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_hit_imf = new TH2F("ptm_rap_hit_imf","MCTrack-hit imf;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);

  fa_v1_rap_hit_pip = new TH2F("v1_rap_hit_pip","MCTrack-hit pi-plus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_pim = new TH2F("v1_rap_hit_pim","MCTrack-hit pi-minus;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_kp  = new TH2F("v1_rap_hit_kp", "MCTrack-hit k-plus;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_km  = new TH2F("v1_rap_hit_km", "MCTrack-hit k-minus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_p   = new TH2F("v1_rap_hit_p",  "MCTrack-hit proton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_pbar= new TH2F("v1_rap_hit_pbar","MCTrack-hit antiproton;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_d   = new TH2F("v1_rap_hit_d",  "MCTrack-hit deuteron;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_t   = new TH2F("v1_rap_hit_t",  "MCTrack-hit triton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_h   = new TH2F("v1_rap_hit_h",  "MCTrack-hit 3he;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_a   = new TH2F("v1_rap_hit_a",  "MCTrack-hit alpha;   y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_hit_imf = new TH2F("v1_rap_hit_imf","MCTrack-hit imf;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  fa_v2_rap_hit_pip = new TH2F("v2_rap_hit_pip","MCTrack-hit pi-plus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_pim = new TH2F("v2_rap_hit_pim","MCTrack-hit pi-minus;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_kp  = new TH2F("v2_rap_hit_kp", "MCTrack-hit k-plus;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_km  = new TH2F("v2_rap_hit_km", "MCTrack-hit k-minus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_p   = new TH2F("v2_rap_hit_p",  "MCTrack-hit proton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_pbar= new TH2F("v2_rap_hit_pbar","MCTrack-hit antiproton;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_d   = new TH2F("v2_rap_hit_d",  "MCTrack-hit deuteron;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_t   = new TH2F("v2_rap_hit_t",  "MCTrack-hit triton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_h   = new TH2F("v2_rap_hit_h",  "MCTrack-hit 3he;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_a   = new TH2F("v2_rap_hit_a",  "MCTrack-hit alpha;   y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_hit_imf = new TH2F("v2_rap_hit_imf","MCTrack-hit imf;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  // GlobalTrack level

  fa_ptm_rap_sts_pip = new TH2F("ptm_rap_sts_pip","MCTrack-sts pi-plus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_pim = new TH2F("ptm_rap_sts_pim","MCTrack-sts pi-minus;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_kp  = new TH2F("ptm_rap_sts_kp", "MCTrack-sts k-plus;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_km  = new TH2F("ptm_rap_sts_km", "MCTrack-sts k-minus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_p   = new TH2F("ptm_rap_sts_p",  "MCTrack-sts proton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_pbar= new TH2F("ptm_rap_sts_pbar","MCTrack-sts antiproton;y;p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_d   = new TH2F("ptm_rap_sts_d",  "MCTrack-sts deuteron;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_t   = new TH2F("ptm_rap_sts_t",  "MCTrack-sts triton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_h   = new TH2F("ptm_rap_sts_h",  "MCTrack-sts 3he;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_a   = new TH2F("ptm_rap_sts_a",  "MCTrack-sts alpha;   y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_sts_imf = new TH2F("ptm_rap_sts_imf","MCTrack-sts imf;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);

  fa_ptm_rap_glo_pip = new TH2F("ptm_rap_glo_pip","MCTrack-glo pi-plus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_pim = new TH2F("ptm_rap_glo_pim","MCTrack-glo pi-minus;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_kp  = new TH2F("ptm_rap_glo_kp", "MCTrack-glo k-plus;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_km  = new TH2F("ptm_rap_glo_km", "MCTrack-glo k-minus; y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_p   = new TH2F("ptm_rap_glo_p",  "MCTrack-glo proton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_pbar= new TH2F("ptm_rap_glo_pbar","MCTrack-glo antiproton;y;p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_d   = new TH2F("ptm_rap_glo_d",  "MCTrack-glo deuteron;y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_t   = new TH2F("ptm_rap_glo_t",  "MCTrack-glo triton;  y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_h   = new TH2F("ptm_rap_glo_h",  "MCTrack-glo 3he;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_a   = new TH2F("ptm_rap_glo_a",  "MCTrack-glo alpha;   y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);
  fa_ptm_rap_glo_imf = new TH2F("ptm_rap_glo_imf","MCTrack-glo imf;     y; p_{T}/m",ptm_nbx,ymin,ymax,ptm_nby,0.,ptmmax);

  fa_v1_rap_glo_pip = new TH2F("v1_rap_glo_pip","MCTrack-glo pi-plus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_pim = new TH2F("v1_rap_glo_pim","MCTrack-glo pi-minus;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_kp  = new TH2F("v1_rap_glo_kp", "MCTrack-glo k-plus;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_km  = new TH2F("v1_rap_glo_km", "MCTrack-glo k-minus; y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_p   = new TH2F("v1_rap_glo_p",  "MCTrack-glo proton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_pbar= new TH2F("v1_rap_glo_pbar","MCTrack-glo antiproton;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_d   = new TH2F("v1_rap_glo_d",  "MCTrack-glo deuteron;y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_t   = new TH2F("v1_rap_glo_t",  "MCTrack-glo triton;  y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_h   = new TH2F("v1_rap_glo_h",  "MCTrack-glo 3he;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_a   = new TH2F("v1_rap_glo_a",  "MCTrack-glo alpha;   y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v1_rap_glo_imf = new TH2F("v1_rap_glo_imf","MCTrack-glo imf;     y; v_{1}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  fa_v2_rap_glo_pip = new TH2F("v2_rap_glo_pip","MCTrack-glo pi-plus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_pim = new TH2F("v2_rap_glo_pim","MCTrack-glo pi-minus;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_kp  = new TH2F("v2_rap_glo_kp", "MCTrack-glo k-plus;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_km  = new TH2F("v2_rap_glo_km", "MCTrack-glo k-minus; y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_p   = new TH2F("v2_rap_glo_p",  "MCTrack-glo proton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_pbar= new TH2F("v2_rap_glo_pbar","MCTrack-glo antiproton;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_d   = new TH2F("v2_rap_glo_d",  "MCTrack-glo deuteron;y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_t   = new TH2F("v2_rap_glo_t",  "MCTrack-glo triton;  y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_h   = new TH2F("v2_rap_glo_h",  "MCTrack-glo 3he;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_a   = new TH2F("v2_rap_glo_a",  "MCTrack-glo alpha;   y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);
  fa_v2_rap_glo_imf = new TH2F("v2_rap_glo_imf","MCTrack-glo imf;     y; v_{2}",v1_nbx,-yvmax,yvmax,v1_nby,-1.,1.);

  // xy - hit densities and rates
  Int_t nbinx=1500;
  Int_t nbiny=1000;
  Float_t xrange=750.;
  Float_t yrange=500.;
  fwxy2 = nbinx*nbiny/4./xrange/yrange;

  fa_xy_poi1 = new TH2F("xy_poi1","TofPoint; ;",nbinx,-xrange,xrange,nbiny,-yrange,yrange);
  fa_xy_poi2 = new TH2F("xy_poi2","TofPoint /cm^{2}; ;",nbinx,-xrange,xrange,nbiny,-yrange,yrange);
  fa_xy_poi3 = new TH2F("xy_poi3","TofPoint /cm^{2}/s; ;",nbinx,-xrange,xrange,nbiny,-yrange,yrange);

  fa_pv_poi = new TH2F("pv_poi","TofPoint(all); velocity;momentum;",100,0.,32.,100.,0.,5.);
  fa_tm_poi = new TH2F("tm_poi","Tofpoi(all); momentum; Tofmass;",100,0.,10.,100.,0.,2.5);
  fa_tm_poiprim = new TH2F("tm_poiprim","Tofpoi(primary); momentum; Tofmass",100,0.,10.,100.,0.,2.5);

  fa_xy_hit1 = new TH2F("xy_hit1","TofHit; ;",nbinx,-xrange,xrange,nbiny,-yrange,yrange);
  fa_xy_hit2 = new TH2F("xy_hit2","TofHit /cm^{2}; ;",nbinx,-xrange,xrange,nbiny,-yrange,yrange);
  fa_xy_hit3 = new TH2F("xy_hit3","TofHit /cm^{2}/s; ;",nbinx,-xrange,xrange,nbiny,-yrange,yrange);


  fa_tof_hit = new TH1F("tof_hit","TofHit(all); t (ns); counts;",100,10.,110.);
  fa_tof_hitprim = new TH1F("tof_hitprim","TofHit(prim); t (ns); counts;",100,10.,110.);

     // PAL Addenda
  fa_tof_hit_res          = new TH1F("tof_hitres",        "TofHit(all); tPt - tHit (ns); counts;",10000,-5.,5.);
  fa_tof_hitprim_res      = new TH1F("tof_hitresprim",    "TofHit(primary); tPt - tHit (ns); counts;", 10000,-5.,5.);
  fa_tof_hit_res_sing     = new TH1F("tof_hitressing",    "TofHit(Single MC); tPt - tHit (ns); counts;",10000,-5.,5.);
  fa_tof_hitprim_res_sing = new TH1F("tof_hitresprimsing","TofHit(primary+singleMC); tPt - tHit (ns); counts;",10000,-5.,5.);
  fa_tof_missmatch_flag   = new TH1F("tof_missmatch_flag","TofHit flag in case of MC Track missmatch; flag []; counts;",100,0.5,100.5);
  fa_tof_missmatch_typ    = new TH1F("tof_missmatch_typ", "Track type in case of MC Track missmatch; ; counts;", 2,-0.5,1.5);
  fa_tof_missmatch_typ->GetXaxis()->SetBinLabel(1,  "Primary");
  fa_tof_missmatch_typ->GetXaxis()->SetBinLabel(2,  "Secondary");
  fa_tof_missmatch_pos    = new TH2F("tof_missmatch_pos", "TofHit pos. in case of MC Track missmatch; X [cm]; Y [cm];",
        nbinx,-xrange,xrange,nbiny,-yrange,yrange);

  fa_pv_hit = new TH2F("pv_hit","TofHit(all); velocity; momentum;",100,0.,32.,100.,0.,5.);
  fa_tm_hit = new TH2F("tm_hit","TofHit(all); momentum; Tofmass;",100,0.,10.,100.,0.,2.5);
  fa_tm_hitprim = new TH2F("tm_hitprim","TofHit(primary); momentum; Tofmass",100,0.,10.,100.,0.,2.5);

  fa_dxx    = new TH2F("dxx","TofHit; x; Delta x;",100,-xrange,xrange,500.,-10.,10.);
  fa_dxy    = new TH2F("dxy","TofHit; y; Delta x;",100,-yrange,yrange,500.,-10.,10.);
//  fa_dxz    = new TH2F("dxz","TofHit; z; Delta x;",100,400.,650.,50.,-10.,10.);
  fa_dxz    = new TH2F("dxz","TofHit; z; Delta x;",100,950.,1200.,500.,-10.,10.);
  fa_dyy    = new TH2F("dyy","TofHit; y; Delta y;",100,-yrange,yrange,500.,-10.,10.);
  fa_dzx    = new TH2F("dzx","TofHit; x; Delta z;",100,-xrange,xrange,500,-20.,20.);
  fa_dzy    = new TH2F("dzy","TofHit; y; Delta z;",100,-yrange,yrange,500,-20.,20.);
//  fa_dzz    = new TH2F("dzz","TofHit; z; Delta z;",100,400.,650.,50,-20.,20.);
  fa_dzz    = new TH2F("dzz","TofHit; z; Delta z;",100,950.,1200.,500,-200.,20.);

  fa_hit_ch = new TH1F("hit_ch","TofHits; channel; rate (Hz/s);",50000,0.,50000.);
  fa_dhit_ch= new TH2F("dhit_ch","Tof Double Hits; channel; counts;",50000,0.,50000.,2,0.5,2.5);

  fa_xy_glo1 = new TH2F("xy_glo1","GlobalTrack(all); x (cm); y(cm);",nbinx,-xrange,xrange,nbiny,-yrange,yrange);

  fa_pv_glo = new TH2F("pv_glo","GlobalTrack(all); velocity; momentum;",100,0.,32.,100.,0.,5.);
  fa_tm_glo = new TH2F("tm_glo","GlobalTrack(all); momentum; Tofmass;",100,0.,10.,100.,0.,2.5);
  fa_chi2_mom_glo = new TH2F("chi2_mom_glo","GlobalTrack(all); momentum; chi2;",100,0.,10.,100.,0.,25.);
  fa_len_mom_glo  = new TH2F("len_mom_glo","GlobalTrack(all); momentum; len;",100,0.,10.,100.,0.,1500.);
  fa_tm_gloprim = new TH2F("tm_gloprim","GlobalTrack(primary); momentum; Tofmass",100,0.,10.,100.,0.,2.5);

  //centrality
  fa_mul_b_gen = new TH2F("mul_b_gen","Centrality - gen;impact parameter b; multiplicity;",15,0.,15.,100,0.,10000.);
  fa_mul_b_poi = new TH2F("mul_b_poi","Centrality - poi;impact parameter b; multiplicity;",15,0.,15.,100,0.,10000.);
  fa_mul_b_hit = new TH2F("mul_b_hit","Centrality - hit;impact parameter b; multiplicity;",15,0.,15.,100,0.,2000.);
  fa_mul_b_glo = new TH2F("mul_b_glo","Centrality - glo;impact parameter b; multiplicity;",15,0.,15.,100,0.,2000.);
  fa_mul_b_had = new TH2F("mul_b_had","Centrality - had;impact parameter b; multiplicity;",15,0.,15.,100,0.,2000.);

  //reaction plane
  fa_cdrp_b_gen = new TH2F("cdrp_b_gen","reaction plane resolution - gen;impact parameter b; cos#Delta#phi_{rp};",15,0.,15.,20,-1.,1.);
  fa_drp_b_gen  = new TH2F("drp_b_gen","#Delta#phi-reaction plane - gen ;impact parameter b;",15,0.,15.,36,-180.,180.);
  fa_phirp_b_gen = new TH2F("phirp_b_gen","#phi_{reaction plane} - gen   ;impact parameter b;",15,0.,15.,36,-180.,180.);
  fa_phgrp_b_gen = new TH2F("phgrp_b_gen","#phi_{G reaction plane} - gen ;impact parameter b;",15,0.,15.,36,-180.,180.);
  fa_phphrp_gen = new TH2F("phphrp_gen","#phi_#phi - gen ;#phi_{rec}; #phi_{gen} ;",36,-180.,180.,36,-180.,180.);
  fa_delrp_b_gen = new TH2F("delrp_b_gen","#Delta#phi_{G}-reaction plane - gen ;impact parameter b;#phi_{rec}-#phi_{gen}",15,0.,15.,36,-180.,180.);
  fa_delrp_b_poi = new TH2F("delrp_b_poi","#Delta#phi_{G}-reaction plane - poi ;impact parameter b;#phi_{rec}-#phi_{gen}",15,0.,15.,36,-180.,180.);
  fa_delrp_b_hit = new TH2F("delrp_b_hit","#Delta#phi_{G}-reaction plane - hit ;impact parameter b;#phi_{rec}-#phi_{gen}",15,0.,15.,36,-180.,180.);
  fa_delrp_b_glo = new TH2F("delrp_b_glo","#Delta#phi_{G}-reaction plane - glo ;impact parameter b;#phi_{rec}-#phi_{gen}",15,0.,15.,36,-180.,180.);

  fa_cdelrp_b_gen = new TH2F("cdelrp_b_gen","reaction plane resolution - gen;impact parameter b;cos(#phi_{rec}-#phi_{gen})",15,0.,15.,20,-1.,1.);
  fa_cdelrp_b_poi = new TH2F("cdelrp_b_poi","reaction plane resolution - poi;impact parameter b;cos(#phi_{rec}-#phi_{gen})",15,0.,15.,20,-1.,1.);
  fa_cdelrp_b_hit = new TH2F("cdelrp_b_hit","reaction plane resolution - hit;impact parameter b;cos(#phi_{rec}-#phi_{gen})",15,0.,15.,20,-1.,1.);
  fa_cdelrp_b_glo = new TH2F("cdelrp_b_glo","reaction plane resolution - glo;impact parameter b;cos(#phi_{rec}-#phi_{gen})",15,0.,15.,20,-1.,1.);
  fa_cdelrp_b_had = new TH2F("cdelrp_b_had","reaction plane resolution - had;impact parameter b;cos(#phi_{rec}-#phi_{gen})",15,0.,15.,20,-1.,1.);

  fa_cdrp_b_poi = new TH2F("cdrp_b_poi","reaction plane resolution - poi;impact parameter b; cos#Delta#phi_{rp};",15,0.,15.,20,-1.,1.);
  fa_drp_b_poi  = new TH2F("drp_b_poi","#Delta#phi-reaction plane - poi ;impact parameter b;",15,0.,15.,36,-180.,180.);

  fa_cdrp_b_hit = new TH2F("cdrp_b_hit","reaction plane resolution - hit;impact parameter b; cos#Delta#phi_{rp};",15,0.,15.,20,-1.,1.);
  fa_drp_b_hit  = new TH2F("drp_b_hit","#Delta#phi-reaction plane - hit ;impact parameter b;",15,0.,15.,36,-180.,180.);

  fa_cdrp_b_glo = new TH2F("cdrp_b_glo","reaction plane resolution - glo;impact parameter b; cos#Delta#phi_{rp};",15,0.,15.,20,-1.,1.);
  fa_drp_b_glo  = new TH2F("drp_b_glo","#Delta#phi-reaction plane - glo ;impact parameter b;",15,0.,15.,36,-180.,180.);

  fa_cdrp_b_had = new TH2F("cdrp_b_had","reaction plane resolution - had;impact parameter b; cos#Delta#phi_{rp};",15,0.,15.,20,-1.,1.);
  fa_drp_b_had  = new TH2F("drp_b_had","#Delta#phi-reaction plane - had ;impact parameter b;",15,0.,15.,36,-180.,180.);

  fa_phirp_gen = new TH1F("phirp_gen","#phi_{reaction plane} - gen   ;#phi_{RPgen};",72,-180.,180.);
  fa_phirp_poi = new TH1F("phirp_poi","#phi_{reaction plane} - poi   ;#phi_{RPpoi};",72,-180.,180.);
  fa_phirp_hit = new TH1F("phirp_hit","#phi_{reaction plane} - hit   ;#phi_{RPhit};",72,-180.,180.);
  fa_phirp_glo = new TH1F("phirp_glo","#phi_{reaction plane} - glo   ;#phi_{RPglo};",72,-180.,180.);
  fa_phirp_had = new TH1F("phirp_had","#phi_{reaction plane} - had   ;#phi_{RPhad};",72,-180.,180.);

  fa_phirps_gen = new TH1F("phirps_gen","#phi_{reaction plane} - gen   ;#phi_{sRPgen};",72,-180.,180.);
  fa_phirps_poi = new TH1F("phirps_poi","#phi_{reaction plane} - poi   ;#phi_{sRPpoi};",72,-180.,180.);
  fa_phirps_hit = new TH1F("phirps_hit","#phi_{reaction plane} - hit   ;#phi_{sRPhit};",72,-180.,180.);
  fa_phirps_glo = new TH1F("phirps_glo","#phi_{reaction plane} - glo   ;#phi_{sRPglo};",72,-180.,180.);

  gDirectory->cd( oldir->GetPath() ); // <= To prevent histos from being sucked in by the param file of the TRootManager!

  cout <<"CbmHadronAnalysis::CreateHistogramms: histograms booked in directory "<<gDirectory->GetName()<<endl;
  }
// ------------------------------------------------------------------
InitStatus CbmHadronAnalysis::ReadPdfFile()
{   // Open PDF file and get histogramms for PID
    if(fPdfFileName.IsNull()) return kSUCCESS;
    TFile* pdfFile = new TFile(fPdfFileName);
    if(NULL == pdfFile) {
	cout << "-E- CbmHadronAnalysis::ReadPdfFile : "
	    << "file " << fPdfFileName << " does not exist!" << endl;
        return kSUCCESS;
    }
    return kSUCCESS;
}

// ------------------------------------------------------------------
InitStatus CbmHadronAnalysis::ReadFlowFile()
{   // Open file and get histogramms for RP corrections 
    if(fFlowFileName.IsNull()) return kSUCCESS;
    fflowFile = new TFile(fFlowFileName);
    if(NULL == fflowFile) {
	cout << "-E- CbmHadronAnalysis::ReadFlowFile : "
	    << "file " << fFlowFileName << " does not exist!" << endl;
        return kSUCCESS;
    }
    cout << "-I- CbmHadronAnalysis::ReadFlowFile : RP corrections from " << fFlowFileName << endl;

    return kSUCCESS;
}

// ------------------------------------------------------------------
InitStatus CbmHadronAnalysis::Init()
{
    // Task initialisation
    FairRootManager* rootMgr = FairRootManager::Instance();
    if(NULL == rootMgr) {
	cout << "-E- CbmHadronAnalysis::Init : "
	    << "ROOT manager is not instantiated." << endl;
        return kFATAL;
    }
    fMCEventHeader = (CbmMCEventHeader*) rootMgr->GetObject("MCEventHeader.");
    if(NULL == fMCEventHeader) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no MC Header Info" << endl;
    }

    fMCTracks = (TClonesArray*) rootMgr->GetObject("MCTrack");
    if(NULL == fMCTracks) {
	cout << "-W- CbmHadronAnalysis::Init : "
            << "no MC track array!" << endl;
    }

    fTofPoints = (TClonesArray*) rootMgr->GetObject("TofPoint");
    if(NULL == fTofPoints) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no TOF point array!" << endl;
    }

    fTofHits = (TClonesArray*) rootMgr->GetObject("TofHit");
    if(NULL == fTofHits) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no TOF Hit array!" << endl;
    }

    fGlobalTracks = (TClonesArray*) rootMgr->GetObject("GlobalTrack");
    if(NULL == fGlobalTracks) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no Global Track array!" << endl;
    }

    fHadrons = (TClonesArray*) rootMgr->GetObject("Hadron");
    if(NULL == fHadrons) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no Hadron array!" << endl;
    }    

    fStsTracks = (TClonesArray*) rootMgr->GetObject("StsTrack");
    if(NULL == fStsTracks) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no STS Track array!" << endl;
    }
    fStsHits = (TClonesArray*) rootMgr->GetObject("StsHit");
    if(NULL == fStsHits) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no STS Hit array!" << endl;
    }
    fStsPoints = (TClonesArray*) rootMgr->GetObject("StsPoint");
    if(NULL == fStsPoints) {
	cout << "-W- CbmHadronAnalysis::Init : "
	    << "no STS Point array!" << endl;
    }

    InitStatus status = ReadPdfFile();
    if(kSUCCESS != status) {
        return status;
    }
    InitStatus statf = ReadFlowFile();
    if(kSUCCESS != statf) {
        return statf;
    }

    Float_t pbeam=GetBeamMomentum();
    {
     Float_t ep = TMath::Sqrt(pbeam*pbeam + M2PROT);
     Float_t gp = ep/TMath::Sqrt(M2PROT);
     Float_t bp = TMath::Sqrt(1. - 1./gp/gp);
     Float_t yp = 0.5*TMath::Log((1.+bp)/(1.-bp));
     SetMidY(yp*0.5);
     cout << "-I- CbmHadronAnalysis: Initialize for pbeam = " << pbeam << " rap:" << GetMidY() << endl;
    } 

    cout << "-I- CbmHadronAnalysis::Init : "
         << "initialisation completed." << endl;

    return kSUCCESS;
}
// ------------------------------------------------------------------



// ------------------------------------------------------------------
void CbmHadronAnalysis::Exec(Option_t* option)
{
    // Task execution

    // Declare variables outside the loop
    CbmMCTrack  *MCTrack;
    CbmStsTrack *StsTrack;
    CbmTofPoint *TofPoint;
    CbmTofHit   *TofHit;
    CbmGlobalTrack *GlobTrack;

    Int_t nMCTracks, nTofPoints, nTofHits, nGlobTracks;
    Int_t pdgCode, Np1, Np2;
    Float_t Qx1, Qy1, Qx2, Qy2, phirp1, phirp2, phirp, delrp, rp_weight;
    Float_t RADDEG=57.29577951;
    Float_t p_MC, px_MC, py_MC, pz_MC;
    Float_t mfrag;

    Int_t TrackP[100000];
    Float_t t_hit;
    Bool_t use_pions_for_flow=kTRUE;

    if(GetBSelMax()>0.) {
      if(fMCEventHeader->GetB()>GetBSelMax()){return;};
    }

    if(GetBSelMin()>0.) {
      if(fMCEventHeader->GetB()<GetBSelMin()){return;};
    }

    Float_t yrp_mid=GetMidY();  // midrapidity -> update from simulation!	

    nMCTracks   = fMCTracks->GetEntriesFast();
    nTofPoints  = fTofPoints->GetEntriesFast();
    nTofHits    = fTofHits->GetEntriesFast();
    nGlobTracks = fGlobalTracks->GetEntriesFast();

    if(0){ //nh-debug 
      cout << "<D> HadronAnalysis::Exec "<<nMCTracks
	   << "," << nTofPoints
	   << "," << nTofHits
	   << "," << nGlobTracks
	   << endl;
      cout << "-D- b = "<<fMCEventHeader->GetB()<< ", phi = " << fMCEventHeader->GetPhi() << endl; 
    }
  
    // generator level 
    fa_mul_b_gen->Fill(fMCEventHeader->GetB(),nMCTracks);

    Qx1=0.; Qy1=0.; Np1=0; Qx2=0.; Qy2=0.; Np2=0; 
    for (Int_t k=0; k< nMCTracks; k++) { // inspect MCTracks 
			
	MCTrack = (CbmMCTrack*) fMCTracks->At(k);
			
	if (MCTrack->GetMotherId()!=-1) continue;    // primary particles only

	// Process only hadrons
	pdgCode = MCTrack->GetPdgCode();
	//        cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<endl;

        if (pdgCode<100000000){
 	 if(211 != TMath::Abs(pdgCode) &&   // pions
	    321 != TMath::Abs(pdgCode) &&   // kaons
	    2212 != TMath::Abs(pdgCode))    // protons
	  continue;	
	} else {
	  mfrag=(pdgCode%1000)/10 * .931494028; // ignoring binding energies ...
	                                   // where is the proper mass stored ? 
	}

        Float_t Phip = RADDEG*atan2(MCTrack->GetPy(),MCTrack->GetPx());
        Float_t dphi = Phip - RADDEG*fMCEventHeader->GetPhi();
	if(dphi<-180.) {dphi +=360.;};
	if(dphi> 180.) {dphi -=360.;};
	dphi = dphi/RADDEG;
	rp_weight = 0.;

        //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
	//  " Mass " << MCTrack->GetMass()<<","<<mfrag<<" Y " << MCTrack->GetRapidity() <<
	//  " Pt " << MCTrack->GetPt() <<endl;

	switch(pdgCode) {
	 case 211 : {
	  fa_ptm_rap_gen_pip->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_gen_pip->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_pip->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
  	 };
	 case -211 : {
	  fa_ptm_rap_gen_pim->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_gen_pim->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_pim->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case 321 : {
	  fa_ptm_rap_gen_kp->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_gen_kp->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_kp->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };
	 case -321 : {
	  fa_ptm_rap_gen_km->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_gen_km->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_km->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };
	case 2212 : { //protons
	  fa_ptm_rap_gen_p->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_gen_p->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_p->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
	 case -2212 : {
	  fa_ptm_rap_gen_pbar->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_gen_pbar->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_pbar->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };

 	 case 1000010020 : {  // deuteron
	  fa_ptm_rap_gen_d->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_gen_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000010030 : {  // triton
	  fa_ptm_rap_gen_t->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_gen_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000020030 : {  // 3he
	  fa_ptm_rap_gen_h->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_gen_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000020040 : {  // alpha
	  fa_ptm_rap_gen_a->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_gen_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
	 
	default : { // intermediate mass fragments 
          //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
	  //" Mass " << MCTrack->GetMass()<<","<<mfrag<<" Y " << MCTrack->GetRapidity() <<
	  //" Pt " << MCTrack->GetPt() <<endl;
	  fa_ptm_rap_gen_imf->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_gen_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_gen_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination  (optimistic)
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
          break;
	};
	}
        if (rp_weight != 0.) {
	   if(gRandom->Uniform(1)>0.5)  { //subdivide events into 2 random subevents
	    Np1++;
            Qx1=Qx1+rp_weight*MCTrack->GetPx();
            Qy1=Qy1+rp_weight*MCTrack->GetPy();
	   }else{
	    Np2++;
            Qx2=Qx2+rp_weight*MCTrack->GetPx();
            Qy2=Qy2+rp_weight*MCTrack->GetPy();	   
	   }
	}
    }

    if (Np1>0 && Np2>0){
      phirp1=atan2(Qy1,Qx1);
      phirp2=atan2(Qy2,Qx2);
      if (fflowFile!=NULL) { // subevent RP flattening 
       TH1F *phirp_gen_fpar = (TH1F *)fflowFile->Get("phirps_gen_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_gen_fpar->GetBinContent(j)  *TMath::Cos(i*phirp1) 
	           +phirp_gen_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp1))/i;
       }
       phirp1+=dphir;

       dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_gen_fpar->GetBinContent(j)  *TMath::Cos(i*phirp2) 
	           +phirp_gen_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp2))/i;
       }
       phirp2+=dphir;
      } // subevent RP flattening end 

      delrp=(phirp1-phirp2);
      fa_phirps_gen->Fill(phirp1*RADDEG);          // 1D histo
      fa_phirps_gen->Fill(phirp2*RADDEG);          // 1D histo
      if(0){ //nh-debug 
       cout << "<D-gen> Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
      }
      fa_cdrp_b_gen->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
      delrp=delrp*RADDEG;
      while(delrp<-180.) {delrp+=360.;}
      while(delrp>180.)  {delrp-=360.;}
      fa_drp_b_gen->Fill(fMCEventHeader->GetB(),delrp);
    }
    phirp=RADDEG*atan2(Qy1+Qy2,Qx1+Qx2);  // full reaction plane
    while(phirp<-180.) {phirp+=360.;}
    while(phirp>180.)  {phirp-=360.;}
    if (fflowFile!=NULL) { // RP flattening 
       TH1F *phirp_gen_fpar = (TH1F *)fflowFile->Get("phirp_gen_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
	 //cout << " RP flat par "<< i << ","<<j<< " par " << phirp_gen_fpar->GetBinContent(j) 
	 //     << ","<< phirp_gen_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
         dphir += (( -phirp_gen_fpar->GetBinContent(j)  *TMath::Cos(i*phirp/RADDEG) 
	             +phirp_gen_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp/RADDEG))/i);
       }
       //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl; 

       phirp+=dphir*RADDEG;
       while(phirp<-180.) {phirp+=360.;}
       while(phirp>180.)  {phirp-=360.;}
    } // RP flattening end 
    delrp=phirp - RADDEG*fMCEventHeader->GetPhi();
    while(delrp<-180.) {delrp+=360.;}
    while(delrp> 180.) {delrp-=360.;}

    fa_phirp_gen->Fill(phirp);
    fa_delrp_b_gen->Fill(fMCEventHeader->GetB(),delrp);
    fa_cdelrp_b_gen->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));
    fa_phirp_b_gen->Fill(fMCEventHeader->GetB(),phirp);
    fa_phgrp_b_gen->Fill(fMCEventHeader->GetB(),RADDEG*fMCEventHeader->GetPhi());
    fa_phphrp_gen->Fill(phirp,RADDEG*fMCEventHeader->GetPhi());


// TofPoint level 

    for (Int_t k=0; k< nMCTracks; k++) { // reset track detected flag
      TrackP[k]=0;
    }
    Qx1=0.; Qy1=0.; Np1=0; Qx2=0.; Qy2=0.; Np2=0; 
    fa_mul_b_poi->Fill(fMCEventHeader->GetB(),nTofPoints);

    for (Int_t l =0; l<nTofPoints; l++) {
      TofPoint = (CbmTofPoint*) fTofPoints->At(l);
      Int_t k = TofPoint->GetTrackID();
      MCTrack = (CbmMCTrack*) fMCTracks->At(k);
      pdgCode = MCTrack->GetPdgCode();
      if (pdgCode>100000000){
	  mfrag=(pdgCode%1000)/10 * .931494028; // ignoring binding energies ...
      }
      px_MC = MCTrack->GetPx();
      py_MC = MCTrack->GetPy();
      pz_MC = MCTrack->GetPz();
      p_MC  = sqrt(px_MC*px_MC+py_MC*py_MC+pz_MC*pz_MC);

      if (TrackP[k]==0) { // for efficiency
        TrackP[k]=1;

        fa_xy_poi1	->Fill(TofPoint->GetX(),TofPoint->GetY());
        fa_xy_poi2	->Fill(TofPoint->GetX(),TofPoint->GetY(),fwxy2);

        Float_t vel = TofPoint->GetLength()/TofPoint->GetTime();
        Float_t bet = vel / clight;  
        if (bet > 0.99999) {bet=0.99999;} 
        Float_t tofmass =p_MC/bet*TMath::Sqrt(1.-bet*bet);
 
        fa_pv_poi->Fill(vel,p_MC);
        fa_tm_poi->Fill(p_MC,tofmass);

	if (MCTrack->GetMotherId()!=-1) continue;    // primary particles only
        fa_tm_poiprim->Fill(p_MC,tofmass);

        Float_t Phip = RADDEG*atan2(MCTrack->GetPy(),MCTrack->GetPx());
        Float_t dphi = Phip - RADDEG*fMCEventHeader->GetPhi();
	if(dphi<-180.) {dphi +=360.;};
	if(dphi> 180.) {dphi -=360.;};
	dphi = dphi/RADDEG;
	rp_weight = 0.;

	switch(pdgCode) {
	 case 211 : {
	  fa_ptm_rap_poi_pip->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case -211 : {
	  fa_ptm_rap_poi_pim->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case 321 : {
	  fa_ptm_rap_poi_kp->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case -321 : {
	  fa_ptm_rap_poi_km->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case 2212 : {
	  fa_ptm_rap_poi_p->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = +1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
	 case -2212 : {
	  fa_ptm_rap_poi_pbar->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };

 	 case 1000010020 : {  // deuteron
	  fa_ptm_rap_poi_d->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_poi_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_poi_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000010030 : {  // triton
	  fa_ptm_rap_poi_t->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_poi_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_poi_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000020030 : {  // 3he
	  fa_ptm_rap_poi_h->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_poi_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_poi_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000020040 : {  // alpha
	  fa_ptm_rap_poi_a->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_poi_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_poi_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
 	 default : { // intermediate mass fragments 
          //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
	  //" Mass " << MCTrack->GetMass()<<","<<mfrag<<" Y " << MCTrack->GetRapidity() <<
	  //" Pt " << MCTrack->GetPt() <<endl;
	  fa_ptm_rap_poi_imf->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/mfrag);
	  fa_v1_rap_poi_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_poi_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination  (optimistic)
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
          break;
	 };


	}
        if (rp_weight != 0.) {
          if(gRandom->Uniform(1)>0.5)  { //subdivide events into 2 random subevents
	    Np1++;
            Qx1=Qx1+rp_weight*MCTrack->GetPx();
            Qy1=Qy1+rp_weight*MCTrack->GetPy();
	  }else{
	    Np2++;
            Qx2=Qx2+rp_weight*MCTrack->GetPx();
            Qy2=Qy2+rp_weight*MCTrack->GetPy();	   
	  }
	}
      }
    }
    if (Np1>0 && Np2>0){
      phirp1=atan2(Qy1,Qx1);
      phirp2=atan2(Qy2,Qx2);
      if (fflowFile!=NULL) { // subevent RP flattening 
       TH1F *phirp_poi_fpar = (TH1F *)fflowFile->Get("phirps_poi_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_poi_fpar->GetBinContent(j)  *TMath::Cos(i*phirp1) 
	           +phirp_poi_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp1))/i;
       }
       phirp1+=dphir;

       dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_poi_fpar->GetBinContent(j)  *TMath::Cos(i*phirp2) 
	           +phirp_poi_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp2))/i;
       }
       phirp2+=dphir;
      } // subevent RP flattening end 
      delrp=(phirp1-phirp2);
      fa_phirps_poi->Fill(phirp1*RADDEG);          // 1D histo
      fa_phirps_poi->Fill(phirp2*RADDEG);          // 1D histo
      if(0){ //nh-debug 
       cout << "<D-poi> Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
      }
      fa_cdrp_b_poi->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
      delrp=delrp*RADDEG;
      if(delrp<-180.) delrp+=360.;
      if(delrp>180.)  delrp-=360.;
      fa_drp_b_poi->Fill(fMCEventHeader->GetB(),delrp);
    }
 
    phirp=RADDEG*atan2(Qy1+Qy2,Qx1+Qx2);  // full reaction plane
    while(phirp<-180.) {phirp+=360.;}
    while(phirp>180.)  {phirp-=360.;}
    if (fflowFile!=NULL) { // RP flattening 
       TH1F *phirp_poi_fpar = (TH1F *)fflowFile->Get("phirp_poi_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
	 //cout << " RP flat par "<< i << ","<<j<< " par " << phirp_gen_fpar->GetBinContent(j) 
	 //     << ","<< phirp_gen_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
         dphir += (( -phirp_poi_fpar->GetBinContent(j)  *TMath::Cos(i*phirp/RADDEG) 
	             +phirp_poi_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp/RADDEG))/i);
       }
       //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl; 

       phirp+=dphir*RADDEG;
       while(phirp<-180.) {phirp+=360.;}
       while(phirp>180.)  {phirp-=360.;}
    } // RP flattening end 
    delrp=phirp - RADDEG*fMCEventHeader->GetPhi();
    while(delrp<-180.) {delrp+=360.;}
    while(delrp> 180.) {delrp-=360.;}
    fa_phirp_poi->Fill(phirp);          // 1D histo
    fa_delrp_b_poi->Fill(fMCEventHeader->GetB(),delrp);
    fa_cdelrp_b_poi->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));

// TofHit   level 
    //TofOff=-0.12; // why ?
    for (Int_t k=0; k< nMCTracks; k++) { // reset track detected flag
      TrackP[k]=0;
    }
    Qx1=0.; Qy1=0.; Np1=0; Qx2=0.; Qy2=0.; Np2=0; 
    fa_mul_b_hit->Fill(fMCEventHeader->GetB(),nTofHits);

    for (Int_t j =0; j<nTofHits; j++) {
      //cout << "<D-hit> j= " << j << endl;
      TofHit   = (CbmTofHit*) fTofHits->At(j);
      Int_t l = TofHit->GetRefId();
      //cout << "<D-hit> l= " << l << endl;
      TofPoint = (CbmTofPoint*) fTofPoints->At(l);
      Int_t k = TofPoint->GetTrackID();
      //cout << "<D-hit> k= " << k << endl;
      MCTrack = (CbmMCTrack*) fMCTracks->At(k);			
      if(0){ //nh-debug 
	cout << "<D-hit> " << j << "," << l << "," << k << endl; 
      }
      pdgCode = MCTrack->GetPdgCode();
      px_MC = MCTrack->GetPx();
      py_MC = MCTrack->GetPy();
      pz_MC = MCTrack->GetPz();
      p_MC  = sqrt(px_MC*px_MC+py_MC*py_MC+pz_MC*pz_MC);
      if (k>100000){
	cout << "<E-hit> Too many MCTracks " << k << " from " << nMCTracks << endl;
        continue;
      }
      if (TrackP[k]==0) { // for efficiency
	//     if (1) {
        TrackP[k]++;
        t_hit	= TofHit->GetTime();
	Float_t delta_tof = TofPoint->GetTime() - t_hit;
	Float_t delta_x	= TofPoint->GetX() - TofHit->GetX();
	Float_t delta_y	= TofPoint->GetY() - TofHit->GetY();
	Float_t delta_z	= TofPoint->GetZ() - TofHit->GetZ();

	// PAL Addenda
  fa_tof_hit_res->Fill(delta_tof);
  if (MCTrack->GetMotherId()==-1)
     fa_tof_hitprim_res->Fill(delta_tof);
  if( 1 == TofHit->GetFlag() )
  {
     fa_tof_hit_res_sing->Fill(delta_tof);
     if(MCTrack->GetMotherId()==-1)
        fa_tof_hitprim_res_sing ->Fill(delta_tof);
  }

	fa_dxx ->Fill(TofPoint->GetX(),delta_x);
	fa_dxy ->Fill(TofPoint->GetY(),delta_x);
	fa_dxz ->Fill(TofPoint->GetZ(),delta_x);
	fa_dyy ->Fill(TofPoint->GetY(),delta_y);
	fa_dzx ->Fill(TofPoint->GetX(),delta_z);
	fa_dzy ->Fill(TofPoint->GetY(),delta_z);
	fa_dzz ->Fill(TofPoint->GetZ(),delta_z);

	fa_xy_hit1->Fill(TofHit->GetX(),TofHit->GetY());
	fa_xy_hit2->Fill(TofHit->GetX(),TofHit->GetY(),fwxy2);
        fa_hit_ch->Fill(TofHit->GetCh());
        fa_dhit_ch->Fill(TofHit->GetCh(),TofHit->GetFlag());

        Float_t vel = TofPoint->GetLength()/t_hit;
        Float_t bet = vel / clight;  
        if (bet > 0.99999) {bet=0.99999;} 
        Float_t tofmass =p_MC/bet*TMath::Sqrt(1.-bet*bet);
 
        fa_tof_hit->Fill(t_hit);
        fa_pv_hit->Fill(vel,p_MC);
        fa_tm_hit->Fill(p_MC,tofmass);
	
	if (MCTrack->GetMotherId()!=-1) continue;    // primary particles only
        if (TrackP[k]>1) continue; // for efficiency consider only first hit of track
        fa_tm_hitprim->Fill(p_MC,tofmass);
        fa_tof_hitprim->Fill(t_hit);

        Float_t Phip = RADDEG*atan2(MCTrack->GetPy(),MCTrack->GetPx());
        Float_t dphi = Phip - RADDEG*fMCEventHeader->GetPhi();
	if(dphi<-180.) {dphi +=360.;};
	if(dphi> 180.) {dphi -=360.;};
	dphi = dphi/RADDEG;
	
	rp_weight = 0.;
 
	switch(pdgCode) {
	 case 211 : {
	  fa_ptm_rap_hit_pip->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_pip->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_pip->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case -211 : {
	  fa_ptm_rap_hit_pim->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_pim->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_pim->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case 321 : {
	  fa_ptm_rap_hit_kp->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_kp->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_kp->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };
	 case -321 : {
	  fa_ptm_rap_hit_km->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_km->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_km->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };
	 case 2212 : {
	  fa_ptm_rap_hit_p->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_p->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_p->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = +1.;}
           else {
	    rp_weight = -1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case -2212 : {
	  fa_ptm_rap_hit_pbar->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_pbar->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_pbar->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };

 	 case 1000010020 : {  // deuteron
	  fa_ptm_rap_hit_d->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000010030 : {  // triton
	  fa_ptm_rap_hit_t->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

	 case 1000020030 : {  // 3he
	  fa_ptm_rap_hit_h->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
	 case 1000020040 : {  // alpha
	  fa_ptm_rap_hit_a->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };

 	 default : { // intermediate mass fragments 
          //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
	  //" Mass " << MCTrack->GetMass()<<","<<MCTrack->GetMass()<<" Y " << MCTrack->GetRapidity() <<
	  //" Pt " << MCTrack->GetPt() <<endl;
	  fa_ptm_rap_hit_imf->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_hit_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_hit_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination  (optimistic)
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
          break;
	 };

	}

        if (rp_weight != 0.) {
	  if(gRandom->Uniform(1)>0.5)  { //subdivide events into 2 random subevents
	    Np1++;
            Qx1=Qx1+rp_weight*MCTrack->GetPx();
            Qy1=Qy1+rp_weight*MCTrack->GetPy();
	  }else{
	    Np2++;
            Qx2=Qx2+rp_weight*MCTrack->GetPx();
            Qy2=Qy2+rp_weight*MCTrack->GetPy();	   
	  }
	}
      }else {
	if(0){
	cout << "<W> CHA: >=2. hit from track k ="<<k
	     << " Hit# "<< j
	  //   << " Str(Cell) "<< TofHit->GetCell()    << "," << TofPoint->GetCell()
	  //   << " Module "   << TofHit->GetModule()  << "," << TofPoint->GetModule()
          //   << " SM "       << TofHit->GetSModule() << "," << TofPoint->GetSModule()
          //   << " SMType "   << TofHit->GetSMtype()  << "," << TofPoint->GetSMtype()
             <<endl;
	}
      }
    }
    if (Np1>0 && Np2>0){
      phirp1=atan2(Qy1,Qx1);
      phirp2=atan2(Qy2,Qx2);
      if (fflowFile!=NULL) { // subevent RP flattening 
       TH1F *phirp_hit_fpar = (TH1F *)fflowFile->Get("phirps_hit_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_hit_fpar->GetBinContent(j)  *TMath::Cos(i*phirp1) 
	           +phirp_hit_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp1))/i;
       }
       phirp1+=dphir;

       dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_hit_fpar->GetBinContent(j)  *TMath::Cos(i*phirp2) 
	           +phirp_hit_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp2))/i;
       }
       phirp2+=dphir;
      } // subevent RP flattening end 
      fa_phirps_hit->Fill(phirp1*RADDEG);          // 1D histo
      fa_phirps_hit->Fill(phirp2*RADDEG);          // 1D histo
      delrp=(phirp1-phirp2);
      if(0){ //nh-debug 
       cout << "<D-hit> Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
      }
      fa_cdrp_b_hit->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
      delrp=delrp*RADDEG;
      if(delrp<-180.) delrp+=360.;
      if(delrp>180.)  delrp-=360.;
      fa_drp_b_hit->Fill(fMCEventHeader->GetB(),delrp);
    }
    phirp=RADDEG*atan2(Qy1+Qy2,Qx1+Qx2);  // full reaction plane
    while(phirp<-180.) {phirp+=360.;}
    while(phirp>180.)  {phirp-=360.;}
    if (fflowFile!=NULL) { // RP flattening 
       TH1F *phirp_hit_fpar = (TH1F *)fflowFile->Get("phirp_hit_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
	 //cout << " RP flat par "<< i << ","<<j<< " par " << phirp_gen_fpar->GetBinContent(j) 
	 //     << ","<< phirp_gen_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
         dphir += (( -phirp_hit_fpar->GetBinContent(j)  *TMath::Cos(i*phirp/RADDEG) 
	             +phirp_hit_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp/RADDEG))/i);
       }
       //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl; 

       phirp+=dphir*RADDEG;
       while(phirp<-180.) {phirp+=360.;}
       while(phirp>180.)  {phirp-=360.;}
    } // RP flattening end 
    delrp=phirp - RADDEG*fMCEventHeader->GetPhi();
    while(delrp<-180.) {delrp+=360.;}
    while(delrp> 180.) {delrp-=360.;}
    fa_phirp_hit->Fill(phirp);          // 1D histo
    fa_delrp_b_hit->Fill(fMCEventHeader->GetB(),delrp);
    fa_cdelrp_b_hit->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));

// GlobalTrack level 

    Int_t NGTofTrack=0;
    Qx1=0.; Qy1=0.; Np1=0; Qx2=0.; Qy2=0.; Np2=0; 
    fa_mul_b_glo->Fill(fMCEventHeader->GetB(),nGlobTracks);

    for (Int_t i=0; i<nGlobTracks; i++) {    // inspect global tracks 
      GlobTrack = (CbmGlobalTrack*) fGlobalTracks->At(i);
      Int_t s = GlobTrack->GetStsTrackIndex();
      Int_t j = GlobTrack->GetTofHitIndex();
      if (0){ // nh-debug
       const FairTrackParam *tpard = GlobTrack->GetParamFirst();
       Float_t qpf=  tpard->GetQp();     //       
       tpard = GlobTrack->GetParamLast();
       Float_t qpl=  tpard->GetQp();     //       

	cout << "<D> GlobTrack " << i << " Tof "<< j 
           << " Chi2 "<<GlobTrack->GetChi2() 
	   << " NDF  "<<GlobTrack->GetNDF() 
	   << " qp  " << qpf << ","<< qpl
           << " len " <<GlobTrack->GetLength()<<endl; 
      }
      Int_t smc=-1;
      if (s>-1) {   // STS Track analysis 
       StsTrack = (CbmStsTrack *) fStsTracks->At(s);
       Int_t NStsHits = StsTrack->GetNStsHits();
       Int_t NStsMCt=0;
       for(Int_t ih = 0; ih < NStsHits; ih++) {
        CbmStsHit* hit = (CbmStsHit*) fStsHits->At(StsTrack->GetStsHitIndex(ih));
        Int_t sh = hit->GetRefIndex();
        if(sh>-1) {
         CbmStsPoint* poi =  (CbmStsPoint*) fStsPoints->At(sh);
         if(smc!=poi->GetTrackID()){
          smc = poi->GetTrackID();
	  if(smc>-1){
	   NStsMCt++;
	   //	   cout << " STS Track info: " << s<<","<< " #hits: "<< NStsHits << " Hit #:"<<ih<<" Ref "<< sh 
	   //               <<" MCTrackIndex:" << smc << " NStsMCt "<<NStsMCt<< endl;
	  }
	 }
	}else {
//	  cout << "StsPoint not available" << endl;
	}
       } // loop over STS hits finished 

       // analysis of STS tracks   
       if (-1<smc) {
       MCTrack = (CbmMCTrack*) fMCTracks->At(smc);			
       pdgCode = MCTrack->GetPdgCode();
       px_MC = MCTrack->GetPx();
       py_MC = MCTrack->GetPy();
       pz_MC = MCTrack->GetPz();
       p_MC  = sqrt(px_MC*px_MC+py_MC*py_MC+pz_MC*pz_MC);       
       if (MCTrack->GetMotherId()==-1) { // select primaries 
	if (0 == MCTrack->GetMass()) continue;
	switch(pdgCode) {
	 case 211 : {
	  fa_ptm_rap_sts_pip->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case -211 : {
	  fa_ptm_rap_sts_pim->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case 321 : {
	  fa_ptm_rap_sts_kp->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case -321 : {
	  fa_ptm_rap_sts_km->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case 2212 : {
	  fa_ptm_rap_sts_p->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	 };
	 case -2212 : {
	  fa_ptm_rap_sts_pbar->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
 	 case 1000010020 : {  // deuteron
	  fa_ptm_rap_sts_d->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case 1000010030 : {  // triton
	  fa_ptm_rap_sts_t->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case 1000020030 : {  // 3he
	  fa_ptm_rap_sts_h->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
	 case 1000020040 : {  // alpha
	  fa_ptm_rap_sts_a->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  break;
	 };
 	 default : { // intermediate mass fragments 
          //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
	  //" Mass " << MCTrack->GetMass()<<","<<MCTrack->GetMass()<<" Y " << MCTrack->GetRapidity() <<
	  //" Pt " << MCTrack->GetPt() <<endl;
	  fa_ptm_rap_sts_imf->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
          break;
	 };
	}
       }
       } 
      }
      if (j>-1) {
       NGTofTrack++;
       TofHit   = (CbmTofHit*) fTofHits->At(j);
       Int_t l = TofHit->GetRefId();
       TofPoint = (CbmTofPoint*) fTofPoints->At(l);
       Int_t k = TofPoint->GetTrackID();
       if (k!=smc){
//          cout << " Ana: MCTrack mismatch: "<< k <<" - "<<smc<<endl;
          fa_tof_missmatch_flag->Fill(TofHit->GetFlag() );
          fa_tof_missmatch_typ->Fill( -1 == MCTrack->GetMotherId() ? kTRUE: kFALSE );
          fa_tof_missmatch_pos->Fill( TofHit->GetX(), TofHit->GetY() );
       }
       MCTrack = (CbmMCTrack*) fMCTracks->At(k);			
       pdgCode = MCTrack->GetPdgCode();
       px_MC = MCTrack->GetPx();
       py_MC = MCTrack->GetPy();
       pz_MC = MCTrack->GetPz();
       p_MC  = sqrt(px_MC*px_MC+py_MC*py_MC+pz_MC*pz_MC);

       Double_t len = GlobTrack->GetLength();
       const FairTrackParam *tpar = GlobTrack->GetParamFirst();
       //       FairTrackParam *tpar = GlobTrack->GetParamLast();
       if (tpar->GetQp()==0.) {
	 cout << "Invalid momentum for global track "<< i << " Hit " << j <<endl;
         continue;
       }
       Double_t mom = 1./tpar->GetQp();
       if(mom<0.) mom=-mom;
       Float_t vel=len/TofHit->GetTime();
       Float_t bet = vel / clight;  
       if (bet > 0.99999) {bet=0.99999;} 
       Float_t tofmass =mom/bet*sqrt(1.-bet*bet);
       Double_t chi2=(Double_t)(GlobTrack->GetChi2())/(GlobTrack->GetNDF());

       fa_xy_glo1->Fill(TofHit->GetX(),TofHit->GetY());
       fa_pv_glo->Fill(vel,mom);
       fa_tm_glo->Fill(mom,tofmass);
       fa_chi2_mom_glo->Fill(mom,chi2);
       fa_len_mom_glo->Fill(mom,len);

       if (MCTrack->GetMotherId()==-1) { // select primaries 
        fa_tm_gloprim->Fill(mom,tofmass);
        Float_t Phip = RADDEG*atan2(MCTrack->GetPy(),MCTrack->GetPx());
        Float_t dphi = Phip - RADDEG*fMCEventHeader->GetPhi();
	if(dphi<-180.) {dphi +=360.;};
	if(dphi> 180.) {dphi -=360.;};
	dphi = dphi/RADDEG;
	rp_weight = 0.;

	switch(pdgCode) {
	 case 211 : {
	  fa_ptm_rap_glo_pip->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_pip->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_pip->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case -211 : {
	  fa_ptm_rap_glo_pim->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_pim->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_pim->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

          if(use_pions_for_flow && TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight = -1.;}
           else {
	    rp_weight = +1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case 321 : {
	  fa_ptm_rap_glo_kp->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_kp->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_kp->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };
	 case -321 : {
	  fa_ptm_rap_glo_km->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_km->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_km->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  break;
	 };
	 case 2212 : {
	  fa_ptm_rap_glo_p->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_p->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_p->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));


	  // reaction plane determination 
         if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }else{
	    rp_weight = 0.;
	  }
	  break;
	 };
	 case -2212 : {
	  fa_ptm_rap_glo_pbar->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_pbar->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_pbar->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));

	  break;
	 };
 	 case 1000010020 : {  // deuteron
	  fa_ptm_rap_glo_d->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_d->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
	 case 1000010030 : {  // triton
	  fa_ptm_rap_glo_t->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_t->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
	 case 1000020030 : {  // 3he
	  fa_ptm_rap_glo_h->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_h->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
	 case 1000020040 : {  // alpha
	  fa_ptm_rap_glo_a->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_a->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination 
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
	  break;
	 };
 	 default : { // intermediate mass fragments 
          //cout << " Track k="<<k<<", pdgCode = "<<pdgCode<<
	  //" Mass " << MCTrack->GetMass()<<","<<MCTrack->GetMass()<<" Y " << MCTrack->GetRapidity() <<
	  //" Pt " << MCTrack->GetPt() <<endl;
	  fa_ptm_rap_glo_imf->Fill(MCTrack->GetRapidity(),MCTrack->GetPt()/MCTrack->GetMass());
	  fa_v1_rap_glo_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(dphi));
	  fa_v2_rap_glo_imf->Fill((MCTrack->GetRapidity()-GetMidY())/GetMidY(),TMath::Cos(2*dphi));
	  // reaction plane determination  (optimistic)
          if(TMath::Abs((MCTrack->GetRapidity()-yrp_mid)/yrp_mid)>GetDY()
	     &&MCTrack->GetPt()/MCTrack->GetMass()>GetFlowMinPtm()){
           if (MCTrack->GetRapidity()>yrp_mid){ // set weights for reaction plane 
	    rp_weight =  1.;}
           else {
	    rp_weight = -1.;
           }
	  }
          break;
	 };

	}
        if (rp_weight != 0.) {
	  if(gRandom->Uniform(1)>0.5)  { //subdivide events into 2 random subevents
	    Np1++;
            Qx1=Qx1+rp_weight*MCTrack->GetPx();
            Qy1=Qy1+rp_weight*MCTrack->GetPy();
	  }else{
	    Np2++;
            Qx2=Qx2+rp_weight*MCTrack->GetPx();
            Qy2=Qy2+rp_weight*MCTrack->GetPy();	   
	  }
	}
       } 
      }
    }
    if (Np1>0 && Np2>0){
      phirp1=atan2(Qy1,Qx1);
      phirp2=atan2(Qy2,Qx2);
      if (fflowFile!=NULL) { // subevent RP flattening 
       TH1F *phirp_glo_fpar = (TH1F *)fflowFile->Get("phirps_glo_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_glo_fpar->GetBinContent(j)  *TMath::Cos(i*phirp1) 
	           +phirp_glo_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp1))/i;
       }
       phirp1+=dphir;

       dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
         dphir += (-phirp_glo_fpar->GetBinContent(j)  *TMath::Cos(i*phirp2) 
	           +phirp_glo_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp2))/i;
       }
       phirp2+=dphir;
      } // subevent RP flattening end 
      delrp=(phirp1-phirp2);
      fa_phirps_glo->Fill(phirp1*RADDEG);          // 1D histo
      fa_phirps_glo->Fill(phirp2*RADDEG);          // 1D histo
      // cout << " Impact parameter "<<fMCEventHeader->GetB()<< ", delrp = "<< delrp << endl;
      fa_cdrp_b_glo->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp));
      delrp=delrp*RADDEG;
      if(delrp<-180.) delrp+=360.;
      if(delrp>180.)  delrp-=360.;
      fa_drp_b_glo->Fill(fMCEventHeader->GetB(),delrp);
    }
    phirp=RADDEG*atan2(Qy1+Qy2,Qx1+Qx2);  // full reaction plane
    while(phirp<-180.) {phirp+=360.;}
    while(phirp>180.)  {phirp-=360.;}
    if (fflowFile!=NULL) { // RP flattening 
       TH1F *phirp_glo_fpar = (TH1F *)fflowFile->Get("phirp_glo_fpar");
       Float_t dphir=0.;
       for (int j=0; j<4; j++){
	 Float_t i = (float)(j+1);
	 //cout << " RP flat par "<< i << ","<<j<< " par " << phirp_glo_fpar->GetBinContent(j) 
	 //     << ","<< phirp_glo_fpar->GetBinContent(j+4) << " phirp "<<phirp<<" dphir "<< dphir << endl;
         dphir += (( -phirp_glo_fpar->GetBinContent(j)  *TMath::Cos(i*phirp/RADDEG) 
	             +phirp_glo_fpar->GetBinContent(j+4)*TMath::Sin(i*phirp/RADDEG))/i);
       }
       //cout << " phirp " << phirp << " dphir " << dphir*RADDEG << endl; 

       phirp+=dphir*RADDEG;
       while(phirp<-180.) {phirp+=360.;}
       while(phirp>180.)  {phirp-=360.;}
    } // RP flattening end 
    delrp=phirp - RADDEG*fMCEventHeader->GetPhi();
    while(delrp<-180.) {delrp+=360.;}
    while(delrp> 180.) {delrp-=360.;}
    fa_phirp_glo->Fill(phirp);          // 1D histo
    fa_delrp_b_glo->Fill(fMCEventHeader->GetB(),delrp);
    fa_cdelrp_b_glo->Fill(fMCEventHeader->GetB(),TMath::Cos(delrp/RADDEG));

    fa_mul_b_had->Fill(fMCEventHeader->GetB(),NGTofTrack);

// Hadron level 


    if(0 == (fEvents%1)) {
	cout << "-I- CbmHadronAnalysis::Exec : "
             << "event " << fEvents << " processed." << endl;
    }
    fEvents += 1;
}
// ------------------------------------------------------------------



// ------------------------------------------------------------------
void CbmHadronAnalysis::Finish()
{
 // Normalisations
 cout << "CbmHadronAnalysis::Finish up with " << fEvents << " analyzed events " << endl;
 
 Double_t sfe = 1./fEvents;
 Double_t sfac = 1.E7;

 cout << "<I> Normalisation factors " << sfe << "," << sfac << endl;

 fa_hit_ch->Scale(sfe*sfac);

 fa_xy_poi2->Scale(sfe);
 fa_xy_poi3->Add(fa_xy_poi2,fa_xy_poi2,sfac,0.);

 fa_xy_hit2->Scale(sfe);
 fa_xy_hit3->Add(fa_xy_hit2,fa_xy_hit2,sfac,0.);

 // Finish of the task execution

 // WriteHistogramms();
}
// ------------------------------------------------------------------


// ------------------------------------------------------------------
void CbmHadronAnalysis::WriteHistogramms()
{
 // Write histogramms to the file

   TDirectory * oldir = gDirectory; // <= To prevent histos from being sucked in by the param file of the TRootManager!
   gROOT->cd(); // <= To prevent histos from being sucked in by the param file of the TRootManager !

// TFile *fHist = new TFile("data/auaumbias.hst.root","RECREATE");
 TFile *fHist = new TFile("HadronAna.hst.root","RECREATE");
 {
    // To get the histos list from gROOT, otherwise we get the list from the current file,
    // which is luckily empty as we just recreated it
    gROOT->cd();

  TIter next(gDirectory->GetList());
  fHist->cd(); // Go back in the file for the actual writing
  TH1 *h;
  TObject* obj;
  while(obj= (TObject*)next()){
      if(obj->InheritsFrom(TH1::Class())){
         h = (TH1*)obj;
         cout << "Write histo " << h->GetTitle() << endl;
         h->Write();
      }
  }
 }
 fHist->ls();
 fHist->Close();

 gDirectory->cd( oldir->GetPath() ); // <= To prevent histos from being sucked in by the param file of the TRootManager
}
// ------------------------------------------------------------------

ClassImp(CbmHadronAnalysis);