StFlowMephi/FemtoDstAnalyzer_ReCentering.C

/**
 * \brief Example of how to read a file (list of files) using StFemtoEvent classes
 *
 * RunFemtoDstAnalyzer.C is an example of reading FemtoDst format.
 * One can use either FemtoDst file or a list of femtoDst files (inFile.lis or
 * inFile.list) as an input, and preform physics analysis
 *
 * \author Grigory Nigmatkulov
 * \date May 29, 2018
 */

// This is needed for calling standalone classes
#define _VANILLA_ROOT_

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

// ROOT headers
#include "TROOT.h"
#include "TFile.h"
#include "TChain.h"
#include "TVector2.h"
#include "TVector3.h"
#include "TTree.h"
#include "TSystem.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.h"
#include "TProfile2D.h"
#include "TStopwatch.h"

// FemtoDst headers
#include "StFemtoDstReader.h"
#include "StFemtoDst.h"
#include "StFemtoEvent.h"
#include "StFemtoTrack.h"

// Load libraries (for ROOT_VERSTION_CODE >= 393215)
#if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
R__LOAD_LIBRARY(/mnt/pool/rhic/4/parfenovpeter/STAR/build/libStFemtoDst.so)
#endif

#include "Constants.h"

// inFile - is a name of name.FemtoDst.root file or a name
//          of a name.lis(t) files, that contains a list of
//          name1.FemtoDst.root, name2.FemtoDst.root, ... files



//Used functions (see them below)
Bool_t isGoodEvent(StFemtoEvent *const &event);
Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
Double_t GetWeight(StFemtoTrack *const &track);
TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
Int_t GetVzBin(Double_t vtxZ);

//_________________
void FemtoDstAnalyzer_ReCentering(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
                                  const Char_t *outFile = "./oReCenteringTest.root")
{

  std::cout << "Hi! Lets do some physics, Master!" << std::endl;
  TStopwatch timer;
  timer.Start();

#if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
  gSystem->Load("../libStFemtoDst.so");
#endif

  StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
  femtoReader->Init();

  // This is a way if you want to spead up IO
  std::cout << "Explicit read status for some branches" << std::endl;
  femtoReader->SetStatus("*", 0);
  femtoReader->SetStatus("Event", 1);
  femtoReader->SetStatus("Track", 1);
  std::cout << "Status has been set" << std::endl;

  std::cout << "Now I know what to read, Master!" << std::endl;

  if (!femtoReader->chain())
  {
    std::cout << "No chain has been found." << std::endl;
  }
  Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  std::cout << "eventsInTree: " << eventsInTree << std::endl;
  Long64_t events2read = femtoReader->chain()->GetEntries();

  std::cout << "Number of events to read: " << events2read << std::endl;

  //Manage ReCentering output

  TProfile2D *p_Qx_FULL_EP[fNharmonics][fArraySize];
  TProfile2D *p_Qy_FULL_EP[fNharmonics][fArraySize];  
  TProfile2D *p_Qx_FULL_SP[fNharmonics][fArraySize];
  TProfile2D *p_Qy_FULL_SP[fNharmonics][fArraySize];
  TProfile2D *p_Qx_EAST_EP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_EAST_EP[fNharmonics][fArraySize][NEtaGaps];  
  TProfile2D *p_Qx_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qx_WEST_EP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_WEST_EP[fNharmonics][fArraySize][NEtaGaps];  
  TProfile2D *p_Qx_WEST_SP[fNharmonics][fArraySize][NEtaGaps];
  TProfile2D *p_Qy_WEST_SP[fNharmonics][fArraySize][NEtaGaps];

  for (int iHarm=0; iHarm<fNharmonics; iHarm++)
  {
    for (int iVtx=0; iVtx<fArraySize; iVtx++)
    {
      p_Qx_FULL_EP[iHarm][iVtx] = new TProfile2D(Form("p_Q%ix_FULL_EP%i",iHarm,iVtx),Form("p_Q%ix_FULL_EP%i",iHarm,iVtx),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
      p_Qy_FULL_EP[iHarm][iVtx] = new TProfile2D(Form("p_Q%iy_FULL_EP%i",iHarm,iVtx),Form("p_Q%iy_FULL_EP%i",iHarm,iVtx),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
      p_Qx_FULL_SP[iHarm][iVtx] = new TProfile2D(Form("p_Q%ix_FULL_SP%i",iHarm,iVtx),Form("p_Q%ix_FULL_SP%i",iHarm,iVtx),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
      p_Qy_FULL_SP[iHarm][iVtx] = new TProfile2D(Form("p_Q%iy_FULL_SP%i",iHarm,iVtx),Form("p_Q%iy_FULL_SP%i",iHarm,iVtx),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
      for (int iEtaGap=0; iEtaGap<NEtaGaps; iEtaGap++)
      {  
        p_Qx_EAST_EP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%ix_EAST_EP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%ix_EAST_EP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
        p_Qy_EAST_EP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%iy_EAST_EP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%iy_EAST_EP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
        p_Qx_EAST_SP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%ix_EAST_SP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%ix_EAST_SP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
        p_Qy_EAST_SP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%iy_EAST_SP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%iy_EAST_SP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);

        p_Qx_WEST_EP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%ix_WEST_EP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%ix_WEST_EP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
        p_Qy_WEST_EP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%iy_WEST_EP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%iy_WEST_EP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
        p_Qx_WEST_SP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%ix_WEST_SP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%ix_WEST_SP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
        p_Qy_WEST_SP[iHarm][iVtx][iEtaGap] = new TProfile2D(Form("p_Q%iy_WEST_SP%i_gap%i",iHarm,iVtx,iEtaGap),Form("p_Q%iy_WEST_SP%i_gap%i",iHarm,iVtx,iEtaGap),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
      }
    }
  }


  //Q-vectors
  TVector2 QvTMP;

  Double_t weight;
  Int_t VtxSign;

  Long64_t goodEventCounter = 0;

  // Loop over events
  for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
  {

    if ((iEvent+1) % 1000 == 0){
      std::cout << "Working on event #[" << (iEvent + 1)
                << "/" << events2read << "]" << std::endl;
    }

    Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
    if (!readEvent)
    {
      std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
      break;
    }

    // Retrieve femtoDst
    StFemtoDst *dst = femtoReader->femtoDst();

    // Retrieve event information
    StFemtoEvent *event = dst->event();
    if (!event)
    {
      std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
      break;
    }


    // Event selection
    if ( !isGoodEvent(event) ) continue;
    goodEventCounter++;

    TVector3 pVtx = event->primaryVertex();

    // Track analysis
    Int_t nTracks = dst->numberOfTracks();


    // Vertex sign
    VtxSign = GetVzBin(pVtx.Z());
    if (VtxSign == -1) continue;

    // Track loop
    for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
    {

      // Retrieve i-th femto track
      StFemtoTrack *femtoTrack = dst->track(iTrk);

      if (!femtoTrack)
        continue;

      // Must be a primary track
      if (!femtoTrack->isPrimary())
        continue;

   
      //Track selection
      if (!isGoodTrack(femtoTrack, energy, pVtx)) continue;      

      //Fill recentered Q-vectors for all TPC (FULL)
      weight = GetWeight(femtoTrack);


      // 2-nd harmonics
      QvTMP.Set(0.,0.);
      QvTMP = CalcQvector(femtoTrack,2);

      p_Qx_FULL_EP[0][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X(),weight);
      p_Qy_FULL_EP[0][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y(),weight);

      p_Qx_FULL_SP[0][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X());
      p_Qy_FULL_SP[0][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y());
      
      for (int iGap=0; iGap<NEtaGaps; iGap++)
      {
        //fill Qvectors from EAST arm
        if (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < 0.)//cutEtaGap[iGap])
        {
          p_Qx_EAST_EP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X(),weight);
          p_Qy_EAST_EP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y(),weight);

          p_Qx_EAST_SP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X());
          p_Qy_EAST_SP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y());
        }
        //fill Qvectors from WEST arm
        if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 0.)//cutEtaGap[iGap])
        {
          p_Qx_WEST_EP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X(),weight);
          p_Qy_WEST_EP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y(),weight);

          p_Qx_WEST_SP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X());
          p_Qy_WEST_SP[0][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y());
        }
      }

      // 3-rd harmonics
      QvTMP.Set(0.,0.);
      QvTMP = CalcQvector(femtoTrack,3);

      p_Qx_FULL_EP[1][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X(),weight);
      p_Qy_FULL_EP[1][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y(),weight);

      p_Qx_FULL_SP[1][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X());
      p_Qy_FULL_SP[1][VtxSign]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y());
      
      for (int iGap=0; iGap<NEtaGaps; iGap++)
      {
        //fill Qvectors from EAST arm
        if (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < 0.)//cutEtaGap[iGap])
        {
          p_Qx_EAST_EP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X(),weight);
          p_Qy_EAST_EP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y(),weight);

          p_Qx_EAST_SP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X());
          p_Qy_EAST_SP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y());
        }
        //fill Qvectors from WEST arm
        if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 0.)//cutEtaGap[iGap])
        {
          p_Qx_WEST_EP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X(),weight);
          p_Qy_WEST_EP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y(),weight);

          p_Qx_WEST_SP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.X());
          p_Qy_WEST_SP[1][VtxSign][iGap]->Fill((Double_t) event->runId(),(Double_t) event->cent16(),(Double_t) QvTMP.Y());
        }
      }

      // Check if track has TOF signal
      if (femtoTrack->isTofTrack())
      {

      } //if( isTofTrack() )

    } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)

  } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)

  femtoReader->Finish();

  TFile *output = new TFile(outFile,"recreate");

  for (int iHarm=0; iHarm<fNharmonics; iHarm++)
  {
    for (int iVtx=0; iVtx<fArraySize; iVtx++)
    {
      p_Qx_FULL_EP[iHarm][iVtx]->Write();
      p_Qy_FULL_EP[iHarm][iVtx]->Write();
      p_Qx_FULL_SP[iHarm][iVtx]->Write();
      p_Qy_FULL_SP[iHarm][iVtx]->Write();
      for (int iEtaGap=0; iEtaGap<NEtaGaps; iEtaGap++)
      {  
        p_Qx_EAST_EP[iHarm][iVtx][iEtaGap]->Write();
        p_Qy_EAST_EP[iHarm][iVtx][iEtaGap]->Write();
        p_Qx_EAST_SP[iHarm][iVtx][iEtaGap]->Write();
        p_Qy_EAST_SP[iHarm][iVtx][iEtaGap]->Write();

        p_Qx_WEST_EP[iHarm][iVtx][iEtaGap]->Write();
        p_Qy_WEST_EP[iHarm][iVtx][iEtaGap]->Write();
        p_Qx_WEST_SP[iHarm][iVtx][iEtaGap]->Write();
        p_Qy_WEST_SP[iHarm][iVtx][iEtaGap]->Write();
      }
    }
  }
  output->Close();

  std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
            << std::endl;
  std::cout << "Good events: " << goodEventCounter << std::endl;
  timer.Stop();
  timer.Print();
}

Bool_t isGoodEvent(StFemtoEvent *const &event)
{
  if (!event) return false;
  if (event == nullptr) return false;

  if (event->primaryVertex().Perp() > cutVtxR) return false;
  if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy)) return false;

  if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz) return false;

  return true;
}

Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
{
  if (!track) return false;
  // if (!track->isPrimary()) return false;
  if (track->nHitsFit() < cutNhits) return false;
  if (track->nHitsPoss() <= cutNhitsPoss) return false;
  if ((Double_t) track->nHitsFit()/track->nHitsPoss() < cutNhitsRatio) return false;
  if (TMath::Abs(track->eta()) >= cutEta) return false;
  
  if (track->pt() <= cutPtMin.at(_energy)) return false;
  if (track->pt() > cutPtMax) return false;
  if (track->ptot() > cutPMax) return false;
  if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
  return true;
}

Double_t GetWeight(StFemtoTrack *const &track)
{
  Double_t weight;
  if (track->pt() <= cutPtWeightEP)
  {
    weight = track->pt();
  }
  else
  {
    weight = cutPtWeightEP;
  }
  return weight;
}

TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
{
  TVector2 qv(0.,0.);

  qv.Set(TMath::Cos(_harm*track->phi()),TMath::Sin(_harm*track->phi()));

  return qv;
}

Int_t GetVzBin(Double_t vtxZ)
{
  for (Int_t i=0; i<fArraySize;i++)
  {
    if (vtxZ >= VtxBins[i] && vtxZ < VtxBins[i+1]) return i;
  }
  return -1;
}