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

#include "TMath.h"
#include "TH1F.h"
#include "TH2F.h"
#include "FairRootManager.h"
#include "CbmMCEventHeader.h"
#include "TClonesArray.h"

#include "eventPlane.h"
#include "CbmPsdHit.h"
#include "CbmMCTrack.h"
#include "CbmStsTrack.h"
#include "CbmTrackMatch.h"
#include "TRandom.h"
#include "FairTrackParam.h"
#include "TVector3.h"
#include "CbmKFTrack.h"
#include "CbmKFVertex.h"
#include "CbmL1PFFitter.h"
#include "L1Field.h"

#include <vector>

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


// -----   Default constructor   -------------------------------------------
eventPlane::eventPlane() : FairTask("EventPlane",1), fHeader(NULL), fMCEventData(NULL), flistPSDhit(NULL), fCbmPsdEvent(NULL), flistSTSMCtrack(NULL), flistSTSRECOtrack(NULL), flistSTStrackMATCH(NULL), flistPV(NULL), fCbmStsEvent(NULL)
{
    hEmod = new TH2F ("hEmod"," Edep vs module ", 40, -100., 100., 40, -100., 100.);

    hY = new TH1F("hY","mc/reco. track rapidity (cm)", 400, -2., 2.);
    hP = new TH1F("hP","reco. track momentum", 1000, 0., 10.);
    hIP = new TH1F ("hIP"," reco. track impact parameter ", 1000, 0., 10.);
    hChi2toNDF = new TH1F("hChi2toNDF"," reco. track chi2 to NDF ", 1000, 0., 30.);
        
    fdoFlat_PSD_Qcorr = kFALSE;
    fdoFlat_PSD_Barr = kFALSE;
    
    fdoFlat_STS_Qcorr = kFALSE;
    fdoFlat_STS_Barr = kFALSE;

    fileName_Qcorr = "";
    fileName_Barr = "";
    
    fy_proj_cm = 0.;
    fy_cm = 0.;
    fEn = 0.;
    fcut_delta = 0.4;
    fdoSTSreco = kTRUE;
}

eventPlane::eventPlane(const char* name, Int_t verbose, Double_t En) : FairTask(name, verbose), fHeader(NULL), fMCEventData(NULL), flistPSDhit(NULL), fCbmPsdEvent(NULL), flistSTSMCtrack(NULL), flistSTSRECOtrack(NULL), flistSTStrackMATCH(NULL), flistPV(NULL), fCbmStsEvent(NULL)
{
    hEmod = new TH2F ("hEmod"," Edep vs module ", 40, -100., 100., 40, -100., 100.);
    
    hY = new TH1F("hY","mc/reco. track rapidity (cm)", 400, -2., 2.);
    hP = new TH1F("hP","reco. track momentum", 1000, 0., 10.);
    hIP = new TH1F ("hIP"," reco. track impact parameter ", 1000, 0., 10.);
    hChi2toNDF = new TH1F("hChi2toNDF"," reco. track chi2 to NDF ", 1000, 0, 30);

    fdoFlat_PSD_Qcorr = kFALSE;
    fdoFlat_PSD_Barr = kFALSE;
    
    fdoFlat_STS_Qcorr = kFALSE;
    fdoFlat_STS_Barr = kFALSE;

    fileName_Qcorr = "";
    fileName_Barr = "";

    fy_proj_cm = 0.;
    fy_cm = 0.;
    fEn = En;
    fcut_delta = 0.4;
    fdoSTSreco = kTRUE;
}
// -------------------------------------------------------------------------


// -----   Destructor   ----------------------------------------------------
eventPlane::~eventPlane()
{
    if ( fHeader ) {
	fHeader->Delete();
	delete fHeader;
    }
    
    if ( fMCEventData ) {
	fMCEventData->Delete();
	delete fMCEventData;
    }
  
    if ( flistPSDhit ) {
	flistPSDhit->Delete();
	delete flistPSDhit;
    }

    if ( fCbmPsdEvent ) {
	fCbmPsdEvent->Delete();
	delete fCbmPsdEvent;
    }

    if ( flistSTSMCtrack ) {
	flistSTSMCtrack->Delete();
	delete flistSTSMCtrack;
    }
    
    if ( flistSTSRECOtrack ) {
	flistSTSRECOtrack->Delete();
	delete flistSTSRECOtrack;
    }
    
    if ( flistSTStrackMATCH ) {
	flistSTStrackMATCH->Delete();
	delete flistSTStrackMATCH;
    }
    
    if ( flistPV ) {
	flistPV->Delete();
	delete flistPV;
    }
    
    if ( fCbmStsEvent ) {
	fCbmStsEvent->Delete();
	delete fCbmStsEvent;
    }
}
// -------------------------------------------------------------------------


// -----   Public method Init   --------------------------------------------
InitStatus eventPlane::Init()
{
  fPi = TMath::Pi();
  
  // Get RootManager
  FairRootManager* ioman = FairRootManager::Instance();
  if ( ! ioman ) {
    cout << "-E- eventPlane::Init: "
	 << "RootManager not instantised!" << endl;
    return kFATAL;
  }

  fHeader = (CbmMCEventHeader*) ioman->GetObject("MCEventHeader.");
  if ( ! fHeader ) {
      LOG(FATAL) << "-E- eventPlane::Init: No event header!" << FairLogger::endl;
    return kERROR;
  }

  // =========== Get input array
  
  flistPSDhit = (TClonesArray*) ioman->GetObject("PsdHit");
  if ( ! flistPSDhit ) {
      LOG(FATAL) << "-E- eventPlane::Init: No PSD hits array!" << FairLogger::endl;
    return kERROR;
  }

  flistSTSMCtrack = (TClonesArray*) ioman->GetObject("MCTrack");
  if ( ! flistSTSMCtrack ) {
      LOG(FATAL) << "-E- eventPlane::Init: No MC Tracks array!" << FairLogger::endl;
    return kERROR;
  }

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

  flistSTStrackMATCH = (TClonesArray*) ioman->GetObject("StsTrackMatch");
  if ( ! flistSTStrackMATCH ) {
      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;
  }    
  
  // =========== Create and register output array
  
  fMCEventData = new CbmMCEventData("");
  ioman->Register("McEvent", "MCEvent", fMCEventData, kTRUE);  
  
  fCbmPsdEvent = new CbmPsdEventData("");
  ioman->Register("PsdEvent", "PSDEvent", fCbmPsdEvent, kTRUE);
  
  fCbmStsEvent = new CbmStsEventData("");
  ioman->Register("StsEvent", "STSEvent", fCbmStsEvent, kTRUE);

  for (Int_t i=0; i<100; i++) // to change if different number of PSD modules
  {
      fX_mod[i] = 0.;
      fY_mod[i] = 0.;
      fphi_mod[i] = 0.;
      fR_mod[i] = 0.;
      fedep_mod[i] = 0.;      
  }
  finc_mod = 0;   

  PSDcoordinates();

  STSprojRapidityCM();  

  //------- Event plane flattening (2 methods)

  // ==== Q-recentering method

  for (Int_t k = 0; k < 7; k++)
  {
      PSD_shift_Qx[k] = 0.;
      PSD_shift_Qy[k] = 0.;
      PSD_RMS_Qx[k] = 1.;
      PSD_RMS_Qy[k] = 1.;

      STS_harmo1_shift_Qx_sub1[k] = 0.;
      STS_harmo1_shift_Qy_sub1[k] = 0.;
      STS_harmo1_RMS_Qx_sub1[k] = 1.;
      STS_harmo1_RMS_Qy_sub1[k] = 1.;

      STS_harmo1_shift_Qx_sub2[k] = 0.;
      STS_harmo1_shift_Qy_sub2[k] = 0.;
      STS_harmo1_RMS_Qx_sub2[k] = 1.;
      STS_harmo1_RMS_Qy_sub2[k] = 1.;

      STS_harmo1_shift_Qx_full[k] = 0.;
      STS_harmo1_shift_Qy_full[k] = 0.;
      STS_harmo1_RMS_Qx_full[k] = 1.;
      STS_harmo1_RMS_Qy_full[k] = 1.;
      
      STS_harmo2_shift_Qx_sub1[k] = 0.;
      STS_harmo2_shift_Qy_sub1[k] = 0.;
      STS_harmo2_RMS_Qx_sub1[k] = 1.;
      STS_harmo2_RMS_Qy_sub1[k] = 1.;

      STS_harmo2_shift_Qx_sub2[k] = 0.;
      STS_harmo2_shift_Qy_sub2[k] = 0.;
      STS_harmo2_RMS_Qx_sub2[k] = 1.;
      STS_harmo2_RMS_Qy_sub2[k] = 1.;

      STS_harmo2_shift_Qx_full[k] = 0.;
      STS_harmo2_shift_Qy_full[k] = 0.;
      STS_harmo2_RMS_Qx_full[k] = 1.;
      STS_harmo2_RMS_Qy_full[k] = 1.;
  }

  //init for PSD
  if (fdoFlat_PSD_Qcorr == kTRUE && fileName_Qcorr == "")
  {
      LOG(FATAL) << "-E- eventPlane::Init: No input file for event plane flattening (Q-recentering)" << FairLogger::endl;
      return kERROR;
  }
  if (fdoFlat_PSD_Qcorr == kTRUE) PSDinit_flat_Qcorr();

  //init for STS
  if (fdoFlat_STS_Qcorr == kTRUE && fileName_Qcorr == "")
  {
      LOG(FATAL) << "-E- eventPlane::Init: No input file for event plane flattening (Q-recentering)" << FairLogger::endl;
      return kERROR;
  }
  if (fdoFlat_STS_Qcorr == kTRUE) STSinit_flat_Qcorr();

  cout << "==============" << endl;

  // ==== Barrette method

  for (Int_t k = 0; k < 7; k++)
  {
      for (Int_t l = 0; l < 8; l++)
      {
	  PSD_mean_cosphi[k][l] = 0.;
	  PSD_mean_sinphi[k][l] = 0.;

          STS_harmo1_mean_cosphi_sub1[k][l] = 0.;
	  STS_harmo1_mean_sinphi_sub1[k][l] = 0.;

          STS_harmo1_mean_cosphi_sub2[k][l] = 0.;
	  STS_harmo1_mean_sinphi_sub2[k][l] = 0.;

          STS_harmo1_mean_cosphi_full[k][l] = 0.;
	  STS_harmo1_mean_sinphi_full[k][l] = 0.;
	  
	  STS_harmo2_mean_cosphi_sub1[k][l] = 0.;
	  STS_harmo2_mean_sinphi_sub1[k][l] = 0.;

          STS_harmo2_mean_cosphi_sub2[k][l] = 0.;
	  STS_harmo2_mean_sinphi_sub2[k][l] = 0.;

          STS_harmo2_mean_cosphi_full[k][l] = 0.;
	  STS_harmo2_mean_sinphi_full[k][l] = 0.;
      }
  }

  // init for PSD
  if (fdoFlat_PSD_Barr == kTRUE && fileName_Barr == "")
  {
      LOG(FATAL) << "-E- eventPlane::Init: No input file for PSD event plane flattening (Barrette method)" << FairLogger::endl;
      return kERROR;
  }
  if (fdoFlat_PSD_Barr == kTRUE) PSDinit_flat_Barr();

  // init for STS
  if (fdoFlat_STS_Barr == kTRUE && fileName_Barr == "")
  {
      LOG(FATAL) << "-E- eventPlane::Init: No input file for STS event plane flattening (Barrette method)" << FairLogger::endl;
      return kERROR;
  }
  if (fdoFlat_STS_Barr == kTRUE) STSinit_flat_Barr();

  cout << "-I- eventPlane: Intialisation successfull " << kSUCCESS<< endl;
  return kSUCCESS;

}
// -------------------------------------------------------------------------


// -----   Public method Exec   --------------------------------------------
void eventPlane::Exec(Option_t* opt)
{
    if ( ! fHeader ) Fatal("Exec", "No event header");
    if ( ! flistPSDhit ) Fatal("Exec", "No PsdHit array");
    if ( ! flistSTSMCtrack ) Fatal("Exec", "No MCtrack array");
    
    // TO DO: special function to fill McEventData ? (no incidence on fill - once per evt) ... only 2 lines ...
    
    PSDmodEdep();
    PSDfullEdep();
    PSDtransverseMomMeth();

    if (fdoSTSreco == kTRUE) 
    {
      if ( ! flistSTSRECOtrack ) Fatal("Exec", "No reco. track array");
      if ( ! flistSTStrackMATCH ) Fatal("Exec", "No track match array");
      if ( ! flistPV ) Fatal("Exec", "No reco. PV array");
      
      STSRECOtransverseMomMeth();
    }
    else STSMCtransverseMomMeth();
}

// -----   Public method   --------------------------------------------

void eventPlane::WriteHistogramms()
{
    hEmod->Write();
    hY->Write();
    hP->Write();
    hIP->Write();
    hChi2toNDF->Write();
    
    hEmod->Delete();
    hY->Delete();
    hP->Delete();
    hIP->Delete();
    hChi2toNDF->Delete();
}

void eventPlane::Finish()
{
    // Finish of the task execution
    WriteHistogramms();
}

// ======================= PSD member functions

void eventPlane::PSDcoordinates()
{
    // calculate the phi angle and distance to origin of center of modules (X and Y lab-coordinate taken from txt file)
        
    ifstream fxypos(fileName_PsdModCoord);      // to change: file in input -> later: mod coord. in PsdHit?
    for (Int_t i=1; i<100; i++) // to change if different number of PSD modules: use fxypos.end() -> continue? then inc for the other loops on module?
    {
	if (fxypos.eof()) break;
	
	fxypos >> fX_mod[i] >> fY_mod[i];
	cout << "module " << i << ": (X, Y) = (" << fX_mod[i] << ", " << fY_mod[i] << ")" << endl;
	
	finc_mod++;
    }
    fxypos.close();
    
    finc_mod -= 1;    
    cout << "evenPlane::PSDcoordinates : #modules = "<< finc_mod << endl;

    for (Int_t i=1; i<=finc_mod; i++) // to change if different number of PSD modules
    {
	fphi_mod[i] = TMath::ATan2(fY_mod[i], fX_mod[i]);
	fR_mod[i] = TMath::Sqrt(fX_mod[i] * fX_mod[i] + fY_mod[i] * fY_mod[i]);
    }
}

void eventPlane::PSDmodEdep()
{
    // calculate energy deposit / module
  
    for (Int_t i=0; i<100; i++) { fedep_mod[i] = 0.; }  // to change if different number of PSD modules

    Int_t nPSDhits = flistPSDhit->GetEntriesFast();
    cout << "evenPlane::PSDmodEdep : # PSD hits = " << nPSDhits << endl;

    CbmPsdHit* hit = NULL;
    Int_t mod;

    for (Int_t ihit=0; ihit<nPSDhits; ihit++) 
    {
	hit = (CbmPsdHit*) flistPSDhit->At(ihit);
	if ( ! hit ) continue;
	mod = hit->GetModuleID();
	
	//mod++; // TO CHANGE!!
	
	fedep_mod[mod] += hit->GetEdep();
    }
}

void eventPlane::PSDfullEdep()
{
    // calculate full energy deposit (with and without beam hole)
  
    Double_t edep_tot = 0.;
    for (Int_t i=1; i<=finc_mod; i++)  // to change if different number of PSD modules
    {
	//if (i == 1 || i == 7 || i == 25 || i == 43 || i== 49) continue; // reject modules nbr: 1, 7, 25, 43, 49 (central & corner modules)
	edep_tot += fedep_mod[i];
    }

    Double_t edep_tot_noHole = 0.;
    for (Int_t i=1; i<=finc_mod; i++)   // to change if different number of PSD modules
    {
	//if (i == 1 || i == 7 || i == 43 || i== 49) continue; // reject modules nbr: 1, 7, 43, 49 (corner modules)
	edep_tot_noHole += fedep_mod[i];
    }

    fCbmPsdEvent->SetEdep_full(edep_tot, edep_tot_noHole);
}

void eventPlane::PSDtransverseMomMeth()
{
    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 << "----- imp: " << B << " ------" << endl;
    cout << "----- RP angle: " << RP << " ------" << endl;

    // MC stored in seperate branch "McEvent" in output
    fMCEventData->SetB(B);
    fMCEventData->SetRP(RP);
    
    fCbmPsdEvent->SetEvtID(evtID);               

    Double_t Qx = 0.;
    Double_t Qy = 0.;
    Double_t EP = 0.;
    Double_t deltaEP = 0.;

    for (Int_t i=1; i<=finc_mod; i++) // to change if different number of PSD modules
    {
        //if (i == 1 || i == 7 || i == 25 || i == 43 || i== 49) continue; // reject modules nbr: 1, 7, 25, 43, 49 (central & corner modules)
 	  
	fR_mod[i] = 1.;  // CHANGE
	  
	Qx += TMath::Cos(fphi_mod[i]) * fR_mod[i] * fedep_mod[i];
	Qy += TMath::Sin(fphi_mod[i]) * fR_mod[i] * fedep_mod[i];

        hEmod->Fill(fX_mod[i], fY_mod[i], fedep_mod[i]);
    }

    //---- event plane flattening: Q-recentering
    
    if (fdoFlat_PSD_Qcorr == kTRUE)
    {
        Int_t k;

	if (B < 2) k = 0;
	if (B >= 2 && B < 4) k = 1;
	if (B >= 4 && B < 6) k = 2;
	if (B >= 6 && B < 8) k = 3;
	if (B >= 8 && B < 10) k = 4;
	if (B >= 10 && B < 12) k = 5;
	if (B >= 12) k = 6;

	Qx -= PSD_shift_Qx[k];
	Qx /= PSD_RMS_Qx[k];
        
	Qy -= PSD_shift_Qy[k];
	Qy /= PSD_RMS_Qy[k];

	// later: binning in centrality for event plane flattening let to user choice
        /*
        itcorr = shift_Qx.begin();
	Qx -= *(itcorr + k);

        itcorr = RMS_Qx.begin();
	Qx /= *(itcorr + k);

        itcorr = shift_Qy.begin();
	Qy -= *(itcorr + k);

        itcorr = RMS_Qy.begin();
	Qy /= *(itcorr + k);*/

    }    

    EP = TMath::ATan2( Qy, Qx );

    // ---- event plane flattening: Barrette method
    
    if (fdoFlat_PSD_Barr == kTRUE)
    {
	Int_t k, m;

	if (B < 2) k = 0;
	if (B >= 2 && B < 4) k = 1;
	if (B >= 4 && B < 6) k = 2;
	if (B >= 6 && B < 8) k = 3;
	if (B >= 8 && B < 10) k = 4;
	if (B >= 10 && B < 12) k = 5;
	if (B >= 12) k = 6;

	for (m = 0; m < 8; m++) { EP += (2/(m+1)) * ( - PSD_mean_sinphi[k][m] * TMath::Cos( (m+1) * EP ) + PSD_mean_cosphi[k][m] * TMath::Sin( (m+1) * EP ) ); }
    }

    // diff RECO phi - MC phi in proper range [-pi, pi] -> for drawing purposes      
    deltaEP = EP - RP;   

    if ( deltaEP >   fPi ) deltaEP -= 2 * fPi;
    if ( deltaEP < - fPi ) deltaEP += 2 * fPi;

    fCbmPsdEvent->SetQvect(Qx, Qy);
    fCbmPsdEvent->SetEP(EP);
    fCbmPsdEvent->SetdeltaEP(deltaEP);
}

//---- event plane flattening: Q-recentering, init for PSD

void eventPlane::PSDinit_flat_Qcorr()
{
    cout << "-- correction file for PSD EP (Q-recentering): " << fileName_Qcorr << endl;
    
    // later: break run if no file nor tree
    
    TFile *f = new TFile(fileName_Qcorr,"R");
    if (!f) cout << "-E- eventPlane::PSDinit_flat_Qcorr: pb with the correction file -----" << endl;
    TTree* tree = (TTree* ) f->Get("cbmsim");
    if (!tree) cout << "-E- eventPlane::PSDinit_flat_Qcorr: pb with the correction TTree -----" << endl;

    // later: binning in centrality for event plane flattening let to user choice
    //
    //shift_Qx.reserve(fNbins_corr);
    //RMS_Qx.reserve(fNbins_corr);
    //shift_Qy.reserve(fNbins_corr);
    //RMS_Qy.reserve(fNbins_corr);

    TString cut;
    Int_t k, l;
    //Double_t bstep, l;
    //bstep = 12./(fNbins_corr-1);

    for (k = 0; k < 7; k++)
    {
	l = 2 * k;

	if (k < 6) cut = Form("McEvent.fB >= %.1d && McEvent.fB < %.1d", l, l + 2);
	else cut = Form("McEvent.fB >= %.1d", l);

	cout << "-- correction for: " << cut << endl;

	tree->Draw("PsdEvent.fQx", cut);
	TH1F* h1  = (TH1F*) tree->GetHistogram();
	//shift_Qx.push_back(h1->GetMean()); // later: binning in centrality for event plane flattening let to user choice
	//RMS_Qx.push_back(h1->GetRMS());
        PSD_shift_Qx[k] = h1->GetMean();
	PSD_RMS_Qx[k] = h1->GetRMS();

        tree->Draw("PsdEvent.fQy", cut);
	TH1F* h2  = (TH1F*) tree->GetHistogram();
	//shift_Qy.push_back(h2->GetMean()); // later: binning in centrality for event plane flattening let to user choice
	//RMS_Qy.push_back(h2->GetRMS());
        PSD_shift_Qy[k] = h2->GetMean();
	PSD_RMS_Qy[k] = h2->GetRMS();
    }
    
    delete f;
}

//---- event plane flattening: Barrette method, init for PSD

void eventPlane::PSDinit_flat_Barr()
{
    cout << "================================================" << endl;
    cout << "-- correction file for PSD EP (Barrette method): " << fileName_Barr << endl;
    
    // later: break run if no file nor tree
    
    TFile *f = new TFile(fileName_Barr,"R");
    if (!f) cout << "-E- eventPlane::PSDinit_flat_Barr: pb with the correction file -----" << endl;
    TTree* tree = (TTree* ) f->Get("cbmsim");
    if (!tree) cout << "-E- eventPlane::PSDinit_flat_Barr: pb with the correction TTree -----" << endl;

    TString cut, coeff;
    Int_t k, l, m;

    for (k = 0; k < 7; k++)
    {
	l = 2 * k;

	if (k < 6) cut = Form("McEvent.fB >= %.1d && McEvent.fB < %.1d", l, l + 2);
	else cut = Form("McEvent.fB >= %.1d", l);

	cout << "-- correction for: " << cut << endl;

	for (m = 0; m < 8; m++)
	{
	    coeff = Form("TMath::Cos( %.1d * PsdEvent.fEP_RECO )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h1 = (TH1F*) tree->GetHistogram();
	    PSD_mean_cosphi[k][m] = h1->GetMean();

	    coeff = Form("TMath::Sin( %.1d * PsdEvent.fEP_RECO )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h2 = (TH1F*) tree->GetHistogram();
	    PSD_mean_sinphi[k][m] = h2->GetMean();
	}
    }
    
    delete f;
}

// ======================= STS member functions

void eventPlane::STSMCtransverseMomMeth()
{
    // Declaration variables
    TRandom* rand = new TRandom();
    Int_t x;                      // random procedure
    Double_t w;
    Int_t pdgCode, trackID, type, motherID;
    Double_t p, px, py, pz, pt, phi, charge, mass, energy, y;
    CbmMCTrack* track;

    Int_t mult_v1 = 0;
    Int_t mult_v2 = 0;
    Double_t Qx_harmo1_sub1 = 0.;
    Double_t Qx_harmo1_sub2 = 0.;
    Double_t Qx_harmo1_full = 0.;
    Double_t Qy_harmo1_sub1 = 0.;
    Double_t Qy_harmo1_sub2 = 0.;
    Double_t Qy_harmo1_full = 0.;
    Double_t EP_harmo1_sub1 = 0.;
    Double_t EP_harmo1_sub2 = 0.;
    Double_t EP_harmo1_full = 0.;
    Double_t deltaEP_harmo1_sub1 = 0.;
    Double_t deltaEP_harmo1_sub2 = 0.;
    Double_t deltaEP_harmo1_full = 0.;

    Double_t Qx_harmo2_sub1 = 0.;
    Double_t Qx_harmo2_sub2 = 0.;
    Double_t Qx_harmo2_full = 0.;
    Double_t Qy_harmo2_sub1 = 0.;
    Double_t Qy_harmo2_sub2 = 0.;
    Double_t Qy_harmo2_full = 0.;
    Double_t EP_harmo2_sub1 = 0.;
    Double_t EP_harmo2_sub2 = 0.;
    Double_t EP_harmo2_full = 0.;
    Double_t deltaEP_harmo2_sub1 = 0.;
    Double_t deltaEP_harmo2_sub2 = 0.;
    Double_t deltaEP_harmo2_full = 0.;
    
    Double_t v1 = 0.;  
    Double_t v2 = 0.;
    
    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 << "----- imp: " << B << " ------" << endl;
    //cout << "----- RP angle: " << RP << " ------" << endl;

    fMCEventData->SetB(B);
    fMCEventData->SetRP(RP);
    
    fCbmStsEvent->SetEvtID(evtID);       

    Int_t nSTStracks = flistSTSMCtrack->GetEntriesFast();
    cout << "evenPlane::STSMCtransverseMomMeth: # STS MC tracks = " << nSTStracks << endl;

    for (Int_t itrack=0; itrack<nSTStracks; itrack++)
    {
	track = (CbmMCTrack*) flistSTSMCtrack->At(itrack);
	if ( ! track ) 
        {
            cout << "evenPlane::STSMCtransverseMomMeth: no track pointer!" << endl;
	    continue;
	}

	// acceptance
	if ( track->AccSTS() < 4 ) continue;

	type = track->GetPdgCode();

	if (type < 1000000000)
	{
	    charge = track->GetCharge();
	    y = track->GetRapidity();
	}
	else
	{
	    //pdg = 1000000000 + 10*1000*z + 10*a + i;
	    charge = TMath::Floor( ( type - 1000000000 ) / 10000 );
	    mass = TMath::Floor( ( type - 1000000000 -  10000 * charge ) / 10 );
	    p = track->GetP();
	    energy = TMath::Sqrt(mass*0.93827203*mass*0.93827203 + p*p);
	    pz = track->GetPz();
	    y = 0.5 * TMath::Log( ( energy + pz ) / ( energy - pz ) );
	}

	if (charge == 0.) continue; // safety condition

	// ======= secondaries?
	//motherID = track->GetMotherId();
	//if ( motherID > -1 ) continue;  // no secondary tracks	

	// ======= rapidity cut for harmo 1

	y -= fy_cm;
	y /= fy_proj_cm;

	if ( y < - fcut_delta ) w = -1.;
	if ( y > - fcut_delta && y < fcut_delta ) w = 0.;
	if ( y >   fcut_delta ) w = 1.;

        //cout << "!!! STS Y cut = " << fcut_delta << endl;

	if (w != 0.) hY->Fill(y);
	
	pt = track->GetPt();
	px = track->GetPx();
	py = track->GetPy();
	phi = TMath::ATan2(py, px);	

	// pion et proton anti-correles (WA98 observation! has an effect for harmo1, not harmo2)
	if ( TMath::Abs(type) == 211 ) phi += TMath::Pi();
	if ( phi >  TMath::Pi() ) phi -= 2 * TMath::Pi();
	if ( phi < -TMath::Pi() ) phi += 2 * TMath::Pi();

	// test reso TH 
	v1 += w * TMath::Cos(phi - RP);    // for checking resolution with theoretical expecation i.e. <cos(n(EP - RP))> vs vn x sqrt(N)
	if (w != 0.) mult_v1++;
	
	v2 += TMath::Cos(2*(phi - RP));	
	mult_v2++;
	  
	//pt = 1.; // decomment if no pt weight wanted (not recommended)  	
	
	x = rand->Integer(2); // create sub-event 1 & sub-event 2
	if ( x == 0)
	{
	    Qx_harmo1_sub1 += w * pt * TMath::Cos(phi);
	    Qy_harmo1_sub1 += w * pt * TMath::Sin(phi);
	    
	    Qx_harmo2_sub1 += pt * TMath::Cos(2*phi);
	    Qy_harmo2_sub1 += pt * TMath::Sin(2*phi);
	}
	else
	{
	    Qx_harmo1_sub2 += w * pt * TMath::Cos(phi);
	    Qy_harmo1_sub2 += w * pt * TMath::Sin(phi);
	    
	    Qx_harmo2_sub2 += pt * TMath::Cos(2*phi);
	    Qy_harmo2_sub2 += pt * TMath::Sin(2*phi);
	}

	Qx_harmo1_full += w * pt * TMath::Cos(phi);
	Qy_harmo1_full += w * pt * TMath::Sin(phi);

	Qx_harmo2_full += pt * TMath::Cos(2*phi);
	Qy_harmo2_full += pt * TMath::Sin(2*phi);
    }

    v1 /= mult_v1;
    v2 /= mult_v2;
    fCbmStsEvent->Setv1v2(v1, v2);

    //---- event plane flattening: Q-recentering
    
    if (fdoFlat_STS_Qcorr == kTRUE)
    {
        Int_t k;

	if (B < 2) k = 0;
	if (B >= 2 && B < 4) k = 1;
	if (B >= 4 && B < 6) k = 2;
	if (B >= 6 && B < 8) k = 3;
	if (B >= 8 && B < 10) k = 4;
	if (B >= 10 && B < 12) k = 5;
	if (B >= 12) k = 6;

	// first order event plane
	Qx_harmo1_sub1 -= STS_harmo1_shift_Qx_sub1[k];
	Qx_harmo1_sub1 /= STS_harmo1_RMS_Qx_sub1[k];
        
	Qy_harmo1_sub1 -= STS_harmo1_shift_Qy_sub1[k];
	Qy_harmo1_sub1 /= STS_harmo1_RMS_Qy_sub1[k];

        Qx_harmo1_sub2 -= STS_harmo1_shift_Qx_sub2[k];
	Qx_harmo1_sub2 /= STS_harmo1_RMS_Qx_sub2[k];
        
	Qy_harmo1_sub2 -= STS_harmo1_shift_Qy_sub2[k];
	Qy_harmo1_sub2 /= STS_harmo1_RMS_Qy_sub2[k];

        Qx_harmo1_full -= STS_harmo1_shift_Qx_full[k];
	Qx_harmo1_full /= STS_harmo1_RMS_Qx_full[k];
        
	Qy_harmo1_full -= STS_harmo1_shift_Qy_full[k];
	Qy_harmo1_full /= STS_harmo1_RMS_Qy_full[k];
	
	// second order event plane
	Qx_harmo2_sub1 -= STS_harmo2_shift_Qx_sub1[k];
	Qx_harmo2_sub1 /= STS_harmo2_RMS_Qx_sub1[k];
        
	Qy_harmo2_sub1 -= STS_harmo2_shift_Qy_sub1[k];
	Qy_harmo2_sub1 /= STS_harmo2_RMS_Qy_sub1[k];

        Qx_harmo2_sub2 -= STS_harmo2_shift_Qx_sub2[k];
	Qx_harmo2_sub2 /= STS_harmo2_RMS_Qx_sub2[k];
        
	Qy_harmo2_sub2 -= STS_harmo2_shift_Qy_sub2[k];
	Qy_harmo2_sub2 /= STS_harmo2_RMS_Qy_sub2[k];

        Qx_harmo2_full -= STS_harmo2_shift_Qx_full[k];
	Qx_harmo2_full /= STS_harmo2_RMS_Qx_full[k];
        
	Qy_harmo2_full -= STS_harmo2_shift_Qy_full[k];
	Qy_harmo2_full /= STS_harmo2_RMS_Qy_full[k];
    }

    EP_harmo1_sub1 = TMath::ATan2( Qy_harmo1_sub1, Qx_harmo1_sub1 );
    EP_harmo1_sub2 = TMath::ATan2( Qy_harmo1_sub2, Qx_harmo1_sub2 );
    EP_harmo1_full = TMath::ATan2( Qy_harmo1_full, Qx_harmo1_full );

    EP_harmo2_sub1 = TMath::ATan2( Qy_harmo2_sub1, Qx_harmo2_sub1 )/2.;
    EP_harmo2_sub2 = TMath::ATan2( Qy_harmo2_sub2, Qx_harmo2_sub2 )/2.;
    EP_harmo2_full = TMath::ATan2( Qy_harmo2_full, Qx_harmo2_full )/2.;
    
    //---- event plane flattening: Barrette method
    
    if (fdoFlat_STS_Barr == kTRUE)
    {
	Int_t k, m;

	if (B < 2) k = 0;
	if (B >= 2 && B < 4) k = 1;
	if (B >= 4 && B < 6) k = 2;
	if (B >= 6 && B < 8) k = 3;
	if (B >= 8 && B < 10) k = 4;
	if (B >= 10 && B < 12) k = 5;
	if (B >= 12) k = 6;

	for (m = 0; m < 8; m++)
	{
	    // first order event plane
	    EP_harmo1_sub1 += (2/(m+1)) * ( - STS_harmo1_mean_sinphi_sub1[k][m] * TMath::Cos( (m+1) * EP_harmo1_sub1 ) + STS_harmo1_mean_cosphi_sub1[k][m] * TMath::Sin( (m+1) * EP_harmo1_sub1 ) );
	    EP_harmo1_sub2 += (2/(m+1)) * ( - STS_harmo1_mean_sinphi_sub2[k][m] * TMath::Cos( (m+1) * EP_harmo1_sub2 ) + STS_harmo1_mean_cosphi_sub2[k][m] * TMath::Sin( (m+1) * EP_harmo1_sub2 ) );
            EP_harmo1_full += (2/(m+1)) * ( - STS_harmo1_mean_sinphi_full[k][m] * TMath::Cos( (m+1) * EP_harmo1_full ) + STS_harmo1_mean_cosphi_full[k][m] * TMath::Sin( (m+1) * EP_harmo1_full ) );
	    
	    // second order event plane
	    EP_harmo2_sub1 += (2/(m+1)) * ( - STS_harmo2_mean_sinphi_sub1[k][m] * TMath::Cos( (m+1) * EP_harmo2_sub1 ) + STS_harmo2_mean_cosphi_sub1[k][m] * TMath::Sin( (m+1) * EP_harmo2_sub1 ) );
	    EP_harmo2_sub2 += (2/(m+1)) * ( - STS_harmo2_mean_sinphi_sub2[k][m] * TMath::Cos( (m+1) * EP_harmo2_sub2 ) + STS_harmo2_mean_cosphi_sub2[k][m] * TMath::Sin( (m+1) * EP_harmo2_sub2 ) );
            EP_harmo2_full += (2/(m+1)) * ( - STS_harmo2_mean_sinphi_full[k][m] * TMath::Cos( (m+1) * EP_harmo2_full ) + STS_harmo2_mean_cosphi_full[k][m] * TMath::Sin( (m+1) * EP_harmo2_full ) );
	}
    }

    // diff RECO phi - MC phi in proper range -> for drawing purposes  
    
    // first order event plane (sub-event 1, sub-event 2, full event): proper range [-pi, pi]
    deltaEP_harmo1_sub1 = EP_harmo1_sub1 - RP;  
    deltaEP_harmo1_sub2 = EP_harmo1_sub2 - RP;
    deltaEP_harmo1_full = EP_harmo1_full - RP;

    if ( deltaEP_harmo1_sub1 >   fPi ) deltaEP_harmo1_sub1 -= 2 * fPi;
    if ( deltaEP_harmo1_sub1 < - fPi ) deltaEP_harmo1_sub1 += 2 * fPi;

    if ( deltaEP_harmo1_sub2 >   fPi ) deltaEP_harmo1_sub2 -= 2 * fPi;
    if ( deltaEP_harmo1_sub2 < - fPi ) deltaEP_harmo1_sub2 += 2 * fPi;

    if ( deltaEP_harmo1_full >   fPi ) deltaEP_harmo1_full -= 2 * fPi;
    if ( deltaEP_harmo1_full < - fPi ) deltaEP_harmo1_full += 2 * fPi;

    // second order event plane (sub-event 1, sub-event 2, full event): proper range [-pi/2, pi/2]    
    deltaEP_harmo2_sub1 = EP_harmo2_sub1 - RP;  
    deltaEP_harmo2_sub2 = EP_harmo2_sub2 - RP;
    deltaEP_harmo2_full = EP_harmo2_full - RP;

    if ( deltaEP_harmo2_sub1 >   fPi/2 ) deltaEP_harmo2_sub1 -= fPi;
    if ( deltaEP_harmo2_sub1 < - fPi/2 ) deltaEP_harmo2_sub1 += fPi;

    if ( deltaEP_harmo2_sub2 >   fPi/2 ) deltaEP_harmo2_sub2 -= fPi;
    if ( deltaEP_harmo2_sub2 < - fPi/2 ) deltaEP_harmo2_sub2 += fPi;

    if ( deltaEP_harmo2_full >   fPi/2 ) deltaEP_harmo2_full -= fPi;
    if ( deltaEP_harmo2_full < - fPi/2 ) deltaEP_harmo2_full += fPi;

    // fill sts event object
    
    fCbmStsEvent->SetMult(mult_v2);
    fCbmStsEvent->SetQvect(Qx_harmo1_sub1, Qy_harmo1_sub1, Qx_harmo1_sub2, Qy_harmo1_sub2, Qx_harmo1_full, Qy_harmo1_full, Qx_harmo2_sub1, Qy_harmo2_sub1, Qx_harmo2_sub2, Qy_harmo2_sub2, Qx_harmo2_full, Qy_harmo2_full);
    fCbmStsEvent->SetEP(EP_harmo1_sub1, EP_harmo1_sub2, EP_harmo1_full, EP_harmo2_sub1, EP_harmo2_sub2, EP_harmo2_full);
    fCbmStsEvent->SetdeltaEP(deltaEP_harmo1_sub1, deltaEP_harmo1_sub2, deltaEP_harmo1_full, deltaEP_harmo2_sub1, deltaEP_harmo2_sub2, deltaEP_harmo2_full);
    
    delete rand;
}

void eventPlane::STSprojRapidityCM()
{
    // 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- eventPlane::STSprojRapidityCM: Please specify beam kinetic energy when instantiating eventPlane object (needed for rapidity selection using STS, done in c.m. frame)" << 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;
}

void eventPlane::STSRECOtransverseMomMeth()
{
    // Declaration variables
    TRandom* rand = new TRandom();
    Int_t x;                      // random procedure
    Double_t w;
    Int_t pdgCode, trackID, type, motherID;
    Double_t p, px, py, pz, pt, phi, charge, mass, energy, y;
    Double_t IPx, IPy, IP;
    Double_t chi2;
    Int_t NDF;
    CbmStsTrack* track;
    CbmTrackMatch* match;
    CbmMCTrack* mctrack;
    CbmKFTrack kftrack;

    Int_t mult_v1 = 0;    
    Int_t mult_v2 = 0;
    Double_t Qx_harmo1_sub1 = 0.;
    Double_t Qx_harmo1_sub2 = 0.;
    Double_t Qx_harmo1_full = 0.;
    Double_t Qy_harmo1_sub1 = 0.;
    Double_t Qy_harmo1_sub2 = 0.;
    Double_t Qy_harmo1_full = 0.;
    Double_t EP_harmo1_sub1 = 0.;
    Double_t EP_harmo1_sub2 = 0.;
    Double_t EP_harmo1_full = 0.;
    Double_t deltaEP_harmo1_sub1 = 0.;
    Double_t deltaEP_harmo1_sub2 = 0.;
    Double_t deltaEP_harmo1_full = 0.;

    Double_t Qx_harmo2_sub1 = 0.;
    Double_t Qx_harmo2_sub2 = 0.;
    Double_t Qx_harmo2_full = 0.;
    Double_t Qy_harmo2_sub1 = 0.;
    Double_t Qy_harmo2_sub2 = 0.;
    Double_t Qy_harmo2_full = 0.;
    Double_t EP_harmo2_sub1 = 0.;
    Double_t EP_harmo2_sub2 = 0.;
    Double_t EP_harmo2_full = 0.;
    Double_t deltaEP_harmo2_sub1 = 0.;
    Double_t deltaEP_harmo2_sub2 = 0.;
    Double_t deltaEP_harmo2_full = 0.;
    
    Double_t v1 = 0.;  
    Double_t v2 = 0.;
    
    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 << "----- imp: " << B << " ------" << endl;
    //cout << "----- RP angle: " << RP << " ------" << endl;

    fMCEventData->SetB(B);
    fMCEventData->SetRP(RP);
    
    fCbmStsEvent->SetEvtID(evtID);        
    
    FairTrackParam *trackParam;
    TVector3 momRec;
    TVector3 posRec;
    
    // MC tracks
    Int_t nMCtracks = flistSTSMCtrack->GetEntriesFast();
    cout << "evenPlane::STSRECOtransverseMomMeth: # STS MC tracks = " << nMCtracks << endl;
    
    // RECO tracks
    Int_t nSTStracks = flistSTSRECOtrack->GetEntriesFast();
    cout << "evenPlane::STSRECOtransverseMomMeth: # 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*) flistSTSRECOtrack->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++)
    {
        track = &vRTracks[itrack];
	
	match = (CbmTrackMatch*) flistSTStrackMATCH->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 TOF particle ID
	//
	mctrack = (CbmMCTrack*) flistSTSMCtrack->At(trackID); 
	type = mctrack->GetPdgCode();
	mass = mctrack->GetMass();	
	
	chi2 = track->GetChi2();
	NDF = track->GetNDF();
	chi2 /= NDF;
	//std::cout << "chi2/NDF = " << chi2 << endl;
	hChi2toNDF->Fill(chi2);
	
	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;
	hIP->Fill(IP);
	
	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);	
	//std::cout << "reco. track p = " << p << endl;
	hP->Fill(p);
	
	energy = TMath::Sqrt(mass*0.93827203*mass*0.93827203 + p*p);
	y = 0.5 * TMath::Log( ( energy + pz ) / ( energy - pz ) );

	// ======= rapidity cut for harmo 1

	y -= fy_cm;
	y /= fy_proj_cm;

	if ( y < - fcut_delta ) w = -1.;
	if ( y > - fcut_delta && y < fcut_delta ) w = 0.;
	if ( y >   fcut_delta ) w = 1.;

	if (w != 0.) hY->Fill(y);

        //cout << "!!! STS Y cut = " << fcut_delta << endl;
	
	phi = TMath::ATan2(py, px);	
	
	// pion et proton anti-correles (WA98 observation! has an effect for harmo1, not harmo2)
	if ( TMath::Abs(type) == 211 ) phi += TMath::Pi();
	if ( phi >  TMath::Pi() ) phi -= 2 * TMath::Pi();
	if ( phi < -TMath::Pi() ) phi += 2 * TMath::Pi();

	v1 += w * TMath::Cos(phi - RP);
	if (w != 0.) mult_v1++; 
	v2 += TMath::Cos(2*(phi - RP));	
	mult_v2++;                 
	
	//pt = 1.; // decomment if no pt weight wanted (not recommended)
	
	x = rand->Integer(2);  // create sub-event 1 & sub-event 2
	if ( x == 0)
	{
	    Qx_harmo1_sub1 += w * pt * TMath::Cos(phi);
	    Qy_harmo1_sub1 += w * pt * TMath::Sin(phi);
	    
	    Qx_harmo2_sub1 += pt * TMath::Cos(2*phi);
	    Qy_harmo2_sub1 += pt * TMath::Sin(2*phi);
	}
	else
	{
	    Qx_harmo1_sub2 += w * pt * TMath::Cos(phi);
	    Qy_harmo1_sub2 += w * pt * TMath::Sin(phi);
	    
	    Qx_harmo2_sub2 += pt * TMath::Cos(2*phi);
	    Qy_harmo2_sub2 += pt * TMath::Sin(2*phi);
	}

	Qx_harmo1_full += w * pt * TMath::Cos(phi);
	Qy_harmo1_full += w * pt * TMath::Sin(phi);

	Qx_harmo2_full += pt * TMath::Cos(2*phi);
	Qy_harmo2_full += pt * TMath::Sin(2*phi);
    }
    
    v1 /= mult_v1;
    v2 /= mult_v2;
    fCbmStsEvent->Setv1v2(v1, v2);    
    
    cout << "evenPlane::STSRECOtransverseMomMeth: # STS reco tracks (no ghost) = " << mult_v2 << endl;

    //---- event plane flattening: Q-recentering
    
    if (fdoFlat_STS_Qcorr == kTRUE)
    {
        Int_t k;

	if (B < 2) k = 0;
	if (B >= 2 && B < 4) k = 1;
	if (B >= 4 && B < 6) k = 2;
	if (B >= 6 && B < 8) k = 3;
	if (B >= 8 && B < 10) k = 4;
	if (B >= 10 && B < 12) k = 5;
	if (B >= 12) k = 6;

	// first order event plane
	Qx_harmo1_sub1 -= STS_harmo1_shift_Qx_sub1[k];
	Qx_harmo1_sub1 /= STS_harmo1_RMS_Qx_sub1[k];
        
	Qy_harmo1_sub1 -= STS_harmo1_shift_Qy_sub1[k];
	Qy_harmo1_sub1 /= STS_harmo1_RMS_Qy_sub1[k];

        Qx_harmo1_sub2 -= STS_harmo1_shift_Qx_sub2[k];
	Qx_harmo1_sub2 /= STS_harmo1_RMS_Qx_sub2[k];
        
	Qy_harmo1_sub2 -= STS_harmo1_shift_Qy_sub2[k];
	Qy_harmo1_sub2 /= STS_harmo1_RMS_Qy_sub2[k];

        Qx_harmo1_full -= STS_harmo1_shift_Qx_full[k];
	Qx_harmo1_full /= STS_harmo1_RMS_Qx_full[k];
        
	Qy_harmo1_full -= STS_harmo1_shift_Qy_full[k];
	Qy_harmo1_full /= STS_harmo1_RMS_Qy_full[k];
	
	// second order event plane
	Qx_harmo2_sub1 -= STS_harmo2_shift_Qx_sub1[k];
	Qx_harmo2_sub1 /= STS_harmo2_RMS_Qx_sub1[k];
        
	Qy_harmo2_sub1 -= STS_harmo2_shift_Qy_sub1[k];
	Qy_harmo2_sub1 /= STS_harmo2_RMS_Qy_sub1[k];

        Qx_harmo2_sub2 -= STS_harmo2_shift_Qx_sub2[k];
	Qx_harmo2_sub2 /= STS_harmo2_RMS_Qx_sub2[k];
        
	Qy_harmo2_sub2 -= STS_harmo2_shift_Qy_sub2[k];
	Qy_harmo2_sub2 /= STS_harmo2_RMS_Qy_sub2[k];

        Qx_harmo2_full -= STS_harmo2_shift_Qx_full[k];
	Qx_harmo2_full /= STS_harmo2_RMS_Qx_full[k];
        
	Qy_harmo2_full -= STS_harmo2_shift_Qy_full[k];
	Qy_harmo2_full /= STS_harmo2_RMS_Qy_full[k];
    }

    EP_harmo1_sub1 = TMath::ATan2( Qy_harmo1_sub1, Qx_harmo1_sub1 );
    EP_harmo1_sub2 = TMath::ATan2( Qy_harmo1_sub2, Qx_harmo1_sub2 );
    EP_harmo1_full = TMath::ATan2( Qy_harmo1_full, Qx_harmo1_full );

    EP_harmo2_sub1 = TMath::ATan2( Qy_harmo2_sub1, Qx_harmo2_sub1 )/2.;
    EP_harmo2_sub2 = TMath::ATan2( Qy_harmo2_sub2, Qx_harmo2_sub2 )/2.;
    EP_harmo2_full = TMath::ATan2( Qy_harmo2_full, Qx_harmo2_full )/2.;
    
    //---- event plane flattening: Barrette method
    
    if (fdoFlat_STS_Barr == kTRUE)
    {
	Int_t k, m;

	if (B < 2) k = 0;
	if (B >= 2 && B < 4) k = 1;
	if (B >= 4 && B < 6) k = 2;
	if (B >= 6 && B < 8) k = 3;
	if (B >= 8 && B < 10) k = 4;
	if (B >= 10 && B < 12) k = 5;
	if (B >= 12) k = 6;

	for (m = 0; m < 8; m++)
	{
	    // first order event plane
	    EP_harmo1_sub1 += (2/(m+1)) * ( - STS_harmo1_mean_sinphi_sub1[k][m] * TMath::Cos( (m+1) * EP_harmo1_sub1 ) + STS_harmo1_mean_cosphi_sub1[k][m] * TMath::Sin( (m+1) * EP_harmo1_sub1 ) );
	    EP_harmo1_sub2 += (2/(m+1)) * ( - STS_harmo1_mean_sinphi_sub2[k][m] * TMath::Cos( (m+1) * EP_harmo1_sub2 ) + STS_harmo1_mean_cosphi_sub2[k][m] * TMath::Sin( (m+1) * EP_harmo1_sub2 ) );
            EP_harmo1_full += (2/(m+1)) * ( - STS_harmo1_mean_sinphi_full[k][m] * TMath::Cos( (m+1) * EP_harmo1_full ) + STS_harmo1_mean_cosphi_full[k][m] * TMath::Sin( (m+1) * EP_harmo1_full ) );
	    
	    // second order event plane
	    EP_harmo2_sub1 += (2/(m+1)) * ( - STS_harmo2_mean_sinphi_sub1[k][m] * TMath::Cos( (m+1) * EP_harmo2_sub1 ) + STS_harmo2_mean_cosphi_sub1[k][m] * TMath::Sin( (m+1) * EP_harmo2_sub1 ) );
	    EP_harmo2_sub2 += (2/(m+1)) * ( - STS_harmo2_mean_sinphi_sub2[k][m] * TMath::Cos( (m+1) * EP_harmo2_sub2 ) + STS_harmo2_mean_cosphi_sub2[k][m] * TMath::Sin( (m+1) * EP_harmo2_sub2 ) );
            EP_harmo2_full += (2/(m+1)) * ( - STS_harmo2_mean_sinphi_full[k][m] * TMath::Cos( (m+1) * EP_harmo2_full ) + STS_harmo2_mean_cosphi_full[k][m] * TMath::Sin( (m+1) * EP_harmo2_full ) );
	}
    }

    // diff RECO phi - MC phi in proper range -> for drawing purposes  
    
    // first order event plane (sub-event 1, sub-event 2, full event): proper range [-pi, pi] 
    deltaEP_harmo1_sub1 = EP_harmo1_sub1 - RP;        
    deltaEP_harmo1_sub2 = EP_harmo1_sub2 - RP;
    deltaEP_harmo1_full = EP_harmo1_full - RP;

    if ( deltaEP_harmo1_sub1 >   fPi ) deltaEP_harmo1_sub1 -= 2 * fPi;
    if ( deltaEP_harmo1_sub1 < - fPi ) deltaEP_harmo1_sub1 += 2 * fPi;

    if ( deltaEP_harmo1_sub2 >   fPi ) deltaEP_harmo1_sub2 -= 2 * fPi;
    if ( deltaEP_harmo1_sub2 < - fPi ) deltaEP_harmo1_sub2 += 2 * fPi;

    if ( deltaEP_harmo1_full >   fPi ) deltaEP_harmo1_full -= 2 * fPi;
    if ( deltaEP_harmo1_full < - fPi ) deltaEP_harmo1_full += 2 * fPi;

    // second order event plane (sub-event 1, sub-event 2, full event): proper range [-pi/2, pi/2]
    deltaEP_harmo2_sub1 = EP_harmo2_sub1 - RP;        
    deltaEP_harmo2_sub2 = EP_harmo2_sub2 - RP;
    deltaEP_harmo2_full = EP_harmo2_full - RP;

    if ( deltaEP_harmo2_sub1 >   fPi/2 ) deltaEP_harmo2_sub1 -= fPi;
    if ( deltaEP_harmo2_sub1 < - fPi/2 ) deltaEP_harmo2_sub1 += fPi;

    if ( deltaEP_harmo2_sub2 >   fPi/2 ) deltaEP_harmo2_sub2 -= fPi;
    if ( deltaEP_harmo2_sub2 < - fPi/2 ) deltaEP_harmo2_sub2 += fPi;

    if ( deltaEP_harmo2_full >   fPi/2 ) deltaEP_harmo2_full -= fPi;
    if ( deltaEP_harmo2_full < - fPi/2 ) deltaEP_harmo2_full += fPi;

    // fill sts event object
    fCbmStsEvent->SetMult(mult_v2);
    fCbmStsEvent->SetQvect(Qx_harmo1_sub1, Qy_harmo1_sub1, Qx_harmo1_sub2, Qy_harmo1_sub2, Qx_harmo1_full, Qy_harmo1_full, Qx_harmo2_sub1, Qy_harmo2_sub1, Qx_harmo2_sub2, Qy_harmo2_sub2, Qx_harmo2_full, Qy_harmo2_full);
    fCbmStsEvent->SetEP(EP_harmo1_sub1, EP_harmo1_sub2, EP_harmo1_full, EP_harmo2_sub1, EP_harmo2_sub2, EP_harmo2_full);
    fCbmStsEvent->SetdeltaEP(deltaEP_harmo1_sub1, deltaEP_harmo1_sub2, deltaEP_harmo1_full, deltaEP_harmo2_sub1, deltaEP_harmo2_sub2, deltaEP_harmo2_full);
    
    delete rand;
}

//---- event plane flattening: Q-recentering, init for STS

void eventPlane::STSinit_flat_Qcorr()
{
    cout << "================================================" << endl;
    cout << "-- correction file for STS EP (Q-recentering): " << fileName_Qcorr << endl;
    
    // later: break run if no file nor tree
    
    TFile *f = new TFile(fileName_Qcorr,"R");
    if (!f) cout << "-E- eventPlane::STSinit_flat_Qcorr: pb with the correction file -----" << endl;
    TTree* tree = (TTree* ) f->Get("cbmsim");
    if (!tree) cout << "-E- eventPlane::STSinit_flat_Qcorr: pb with the correction TTree -----" << endl;

    TString cut;
    Int_t k, l;

    for (k = 0; k < 7; k++)
    {
	l = 2 * k;

	if (k < 6) cut = Form("McEvent.fB >= %.1d && McEvent.fB < %.1d", l, l + 2);
	else cut = Form("McEvent.fB >= %.1d", l);

	cout << "-- correction for: " << cut << endl;

	// first order event plane
	tree->Draw("StsEvent.fQx_harmo1_sub1", cut);
	TH1F* h1  = (TH1F*) tree->GetHistogram();
        STS_harmo1_shift_Qx_sub1[k] = h1->GetMean();
	STS_harmo1_RMS_Qx_sub1[k] = h1->GetRMS();

        tree->Draw("StsEvent.fQy_harmo1_sub1", cut);
	TH1F* h2  = (TH1F*) tree->GetHistogram();
        STS_harmo1_shift_Qy_sub1[k] = h2->GetMean();
	STS_harmo1_RMS_Qy_sub1[k] = h2->GetRMS();

        tree->Draw("StsEvent.fQx_harmo1_sub2", cut);
	TH1F* h3  = (TH1F*) tree->GetHistogram();
        STS_harmo1_shift_Qx_sub2[k] = h3->GetMean();
	STS_harmo1_RMS_Qx_sub2[k] = h3->GetRMS();

        tree->Draw("StsEvent.fQy_harmo1_sub2", cut);
	TH1F* h4  = (TH1F*) tree->GetHistogram();
        STS_harmo1_shift_Qy_sub2[k] = h4->GetMean();
	STS_harmo1_RMS_Qy_sub2[k] = h4->GetRMS();

        tree->Draw("StsEvent.fQx_harmo1_full", cut);
	TH1F* h5  = (TH1F*) tree->GetHistogram();
        STS_harmo1_shift_Qx_full[k] = h5->GetMean();
	STS_harmo1_RMS_Qx_full[k] = h5->GetRMS();

        tree->Draw("StsEvent.fQy_harmo1_full", cut);
	TH1F* h6  = (TH1F*) tree->GetHistogram();
        STS_harmo1_shift_Qy_full[k] = h6->GetMean();
	STS_harmo1_RMS_Qy_full[k] = h6->GetRMS();
		
	// second order event plane
	tree->Draw("StsEvent.fQx_harmo2_sub1", cut);
	TH1F* h7  = (TH1F*) tree->GetHistogram();
        STS_harmo2_shift_Qx_sub1[k] = h7->GetMean();
	STS_harmo2_RMS_Qx_sub1[k] = h7->GetRMS();

        tree->Draw("StsEvent.fQy_harmo2_sub1", cut);
	TH1F* h8  = (TH1F*) tree->GetHistogram();
        STS_harmo2_shift_Qy_sub1[k] = h8->GetMean();
	STS_harmo2_RMS_Qy_sub1[k] = h8->GetRMS();

        tree->Draw("StsEvent.fQx_harmo2_sub2", cut);
	TH1F* h9  = (TH1F*) tree->GetHistogram();
        STS_harmo2_shift_Qx_sub2[k] = h9->GetMean();
	STS_harmo2_RMS_Qx_sub2[k] = h9->GetRMS();

        tree->Draw("StsEvent.fQy_harmo2_sub2", cut);
	TH1F* h10  = (TH1F*) tree->GetHistogram();
        STS_harmo2_shift_Qy_sub2[k] = h10->GetMean();
	STS_harmo2_RMS_Qy_sub2[k] = h10->GetRMS();

        tree->Draw("StsEvent.fQx_harmo2_full", cut);
	TH1F* h11  = (TH1F*) tree->GetHistogram();
        STS_harmo2_shift_Qx_full[k] = h11->GetMean();
	STS_harmo2_RMS_Qx_full[k] = h11->GetRMS();

        tree->Draw("StsEvent.fQy_harmo2_full", cut);
	TH1F* h12  = (TH1F*) tree->GetHistogram();
        STS_harmo2_shift_Qy_full[k] = h12->GetMean();
	STS_harmo2_RMS_Qy_full[k] = h12->GetRMS();
    }
    
    delete f;
}

//---- event plane flattening: Barrette method, init for STS

void eventPlane::STSinit_flat_Barr()
{
    cout << "================================================" << endl;
    cout << "-- correction file for STS EP (Barrette method): " << fileName_Barr << endl;
    
    // later: break run if no file nor tree
    
    TFile *f = new TFile(fileName_Barr,"R");
    if (!f) cout << "-E- eventPlane::STSinit_flat_Barr: pb with the correction file -----" << endl;
    TTree* tree = (TTree* ) f->Get("cbmsim");
    if (!tree) cout << "-E- eventPlane::STSinit_flat_Barr: pb with the correction TTree -----" << endl;

    TString cut, coeff;
    Int_t k, l, m;

    for (k = 0; k < 7; k++)
    {
	l = 2 * k;

	if (k < 6) cut = Form("McEvent.fB >= %.1d && McEvent.fB < %.1d", l, l + 2);
	else cut = Form("McEvent.fB >= %.1d", l);

	cout << "-- correction for: " << cut << endl;

	for (m = 0; m < 8; m++)
	{
	    // first order event plane
	    coeff = Form("TMath::Cos( %.1d * StsEvent.fEP_RECO_harmo1_sub1 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h1 = (TH1F*) tree->GetHistogram();
	    STS_harmo1_mean_cosphi_sub1[k][m] = h1->GetMean();

	    coeff = Form("TMath::Sin( %.1d * StsEvent.fEP_RECO_harmo1_sub1 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h2 = (TH1F*) tree->GetHistogram();
	    STS_harmo1_mean_sinphi_sub1[k][m] = h2->GetMean();

            coeff = Form("TMath::Cos( %.1d * StsEvent.fEP_RECO_harmo1_sub2 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h3 = (TH1F*) tree->GetHistogram();
	    STS_harmo1_mean_cosphi_sub2[k][m] = h3->GetMean();

	    coeff = Form("TMath::Sin( %.1d * StsEvent.fEP_RECO_harmo1_sub2 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h4 = (TH1F*) tree->GetHistogram();
	    STS_harmo1_mean_sinphi_sub2[k][m] = h4->GetMean();

            coeff = Form("TMath::Cos( %.1d * StsEvent.fEP_RECO_harmo1_full )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h5 = (TH1F*) tree->GetHistogram();
	    STS_harmo1_mean_cosphi_full[k][m] = h5->GetMean();

	    coeff = Form("TMath::Sin( %.1d * StsEvent.fEP_RECO_harmo1_full )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h6 = (TH1F*) tree->GetHistogram();
	    STS_harmo1_mean_sinphi_full[k][m] = h6->GetMean();
	    	    
	    // second order event plane
	    coeff = Form("TMath::Cos( %.1d * StsEvent.fEP_RECO_harmo2_sub1 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h7 = (TH1F*) tree->GetHistogram();
	    STS_harmo2_mean_cosphi_sub1[k][m] = h7->GetMean();

	    coeff = Form("TMath::Sin( %.1d * StsEvent.fEP_RECO_harmo2_sub1 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h8 = (TH1F*) tree->GetHistogram();
	    STS_harmo2_mean_sinphi_sub1[k][m] = h8->GetMean();

            coeff = Form("TMath::Cos( %.1d * StsEvent.fEP_RECO_harmo2_sub2 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h9 = (TH1F*) tree->GetHistogram();
	    STS_harmo2_mean_cosphi_sub2[k][m] = h9->GetMean();

	    coeff = Form("TMath::Sin( %.1d * StsEvent.fEP_RECO_harmo2_sub2 )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h10 = (TH1F*) tree->GetHistogram();
	    STS_harmo2_mean_sinphi_sub2[k][m] = h10->GetMean();

            coeff = Form("TMath::Cos( %.1d * StsEvent.fEP_RECO_harmo2_full )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h11 = (TH1F*) tree->GetHistogram();
	    STS_harmo2_mean_cosphi_full[k][m] = h11->GetMean();

	    coeff = Form("TMath::Sin( %.1d * StsEvent.fEP_RECO_harmo2_full )", m+1);
	    tree->Draw(coeff, cut);
	    TH1F* h12 = (TH1F*) tree->GetHistogram();
	    STS_harmo2_mean_sinphi_full[k][m] = h12->GetMean();
	}
    }
    
    delete f;
}

ClassImp(eventPlane)