StFemtoEvent/macros/FemtoDstAnalyzer.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>

// ROOT headers
#include "TROOT.h"
#include "TFile.h"
#include "TChain.h"
#include "TTree.h"
#include "TSystem.h"
#include "TH1.h"
#include "TH2.h"
#include "TMath.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(../build/libStFemtoDst.so)
#endif

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

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

  std::cout << "Hi! Lets do some physics, Master!" << std::endl;

  #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;

  // Init output file and histograms
  TFile *fo = new TFile(outFile, "recreate");
  // Vertex position
  TH1D *hVtxX = new TH1D("hVtxX", "hVtxX", 200, -10., 10.);
  TH1D *hVtxY = new TH1D("hVtxY", "hVtxY", 200, -10., 10.);
  TH1D *hVtxZ = new TH1D("hVtxZ", "hVtxZ", 200, -100., 100.);
  TH1D *hVpdZ = new TH1D("hVpdZ", "hVpdZ", 200, -100., 100.);
  // Energy from BBCs
  TH1D *hBbcEastE = new TH1D("hBbcEastE", "hBbcEastE", 200, 0., 2000.);
  TH1D *hBbcWestE = new TH1D("hBbcWestE", "hBbcWestE", 200, 0., 2000.);
  // Energy from ZDC-SMDs
  TH1D *hZdcEastE = new TH1D("hZdcEastE", "hZdcEastE", 200, 0., 2000.);
  TH1D *hZdcWestE = new TH1D("hZdcWestE", "hZdcWestE", 200, 0., 2000.);
  // Centrality (16 bins in total covering 0-80%)
  TH1D *hCent16 = new TH1D("hCent16", "hCent16", 16, 0., 80.);
  // refMults
  TH1D *hRefMult = new TH1D("hRefMult", "hRefMult", 1000., 0., 1000.);
  TH1D *hRefMult2 = new TH1D("hRefMult2", "hRefMult2", 1000., 0., 1000.);
  TH1D *hRefMultCorr = new TH1D("hRefMultCorr", "hRefMultCorr", 1000., 0., 1000.);
  TH1D *hRefMultCorrW = new TH1D("hRefMultCorrW", "hRefMultCorrW", 1000., 0., 1000.);
  TH1D *hNTofHit = new TH1D("hNTofHit", "hNTofHit", 2000., 0., 2000.);
  // Track information
  TH1D *hPt = new TH1D("hPt", "hPt", 300, 0., 3.); // pt
  TH1D *hPhi = new TH1D("hPhi", "hPhi", 360, -1.*TMath::Pi(), TMath::Pi()); // phi
  TH1D *hEta = new TH1D("hEta", "hEta", 240, -1.2, 1.2); // eta
  TH1D *hChi2 = new TH1D("hChi2", "hChi2", 100, 0., 10.); // Chi2 of track (track quality)
  TH1D *hNhits = new TH1D("hNhits", "hNhits", 100, 0., 100.); // N_hits
  TH1D *hNhitsFit = new TH1D("hNhitsFit", "hNhitsFit", 100, 0., 100.); // fitted N_hits
  TH1D *hNhitsPoss = new TH1D("hNhitsPoss", "hNhitsPoss", 100, 0., 100.); // maximum possible N_hits
  TH1D *hDedx = new TH1D("hDedx", "hDedx", 500, 0., 1e-4); // dE/dx
  TH1D *hTofM2 = new TH1D("hTofM2", "hTofM2", 600, -1., 5.); // mass^2 from TOF+TPC
  TH1D *hNsigEl = new TH1D("hNsigEl", "hNsigEl", 200, -10., 10.); // possibility of track being electron calculated from dE/dx in sigmas
  TH1D *hNsigPi = new TH1D("hNsigPi", "hNsigPi", 200, -10., 10.); // possibility of track being pion calculated from dE/dx in sigmas
  TH1D *hNsigKa = new TH1D("hNsigKa", "hNsigKa", 200, -10., 10.); // possibility of track being kaon calculated from dE/dx in sigmas
  TH1D *hNsigPr = new TH1D("hNsigPr", "hNsigPr", 200, -10., 10.); // possibility of track being proton calculated from dE/dx in sigmas
  TH1D *hDCA = new TH1D("hDCA", "hDCA", 100, 0., 5.); // abs value of DCA (sqrt(dcax^2+dcay^2+dcaz^2))
  TH1D *hCharge = new TH1D("hCharge", "hCharge", 4, -2., 2.); // charge

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

    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;
    }

    // Return primary vertex position
    TVector3 pVtx = event->primaryVertex();

    hVtxX->Fill(pVtx.X());
    hVtxY->Fill(pVtx.Y());
    hVtxZ->Fill(pVtx.Z());
    hVpdZ->Fill(event->vpdVz());

    // Reject vertices that are far from the central membrane along the beam
    if( TMath::Abs( pVtx.Z() ) > 200. ) continue;

    // Get centrality from 16 bins:
    /// 15 = 0-5%
    /// 14 = 5-10%
    /// 13 = 10-15%
    /// 12 = 15-20%
    /// 11 = 20-25%
    /// 10 = 25-30%
    ///  9 = 30-35%
    ///  8 = 35-40%
    ///  7 = 40-45%
    ///  6 = 45-50%
    ///  5 = 50-55%
    ///  4 = 55-60%
    ///  3 = 60-65%
    ///  2 = 65-70%
    ///  1 = 70-75%
    ///  0 = 75-80%
    hCent16->Fill( (15.-(double)event->cent16())*5.+2.5 ); // (15 - CentBin)*5 + 2.5 - formula for translating centrality bin to centrality

    hRefMult->Fill(event->refMult());
    hRefMult2->Fill(event->refMult2());
    hRefMultCorr->Fill(event->refMultCorr());
    hRefMultCorrW->Fill(event->refMultCorrWeight());
    hNTofHit->Fill(event->numberOfBTofHit());

    // Fill information from BBCs
    Double_t BbcE_ene = 0.;
    Double_t BbcW_ene = 0.;
    for (int i=0; i<16; i++)
    {
      BbcE_ene += event->bbcAdcEast(i);
      BbcW_ene += event->bbcAdcWest(i);
    }
    hBbcEastE->Fill(BbcE_ene);
    hBbcWestE->Fill(BbcW_ene);

    // Fill information from ZDC-SMDs
    Double_t ZdcE_ene_horizontal = 0.;
    Double_t ZdcE_ene_vertical = 0.;
    Double_t ZdcW_ene_horizontal = 0.;
    Double_t ZdcW_ene_vertical = 0.;
    for (int i=0; i<8; i++)
    {
      ZdcE_ene_horizontal += event->zdcSmdEastHorizontal(i);
      ZdcW_ene_horizontal += event->zdcSmdWestHorizontal(i);
    }
    for (int i=0; i<7; i++)
    {
      ZdcE_ene_vertical += event->zdcSmdEastVertical(i);
      ZdcW_ene_vertical += event->zdcSmdWestVertical(i);
    }
    hZdcEastE->Fill(ZdcE_ene_horizontal+ZdcE_ene_vertical);
    hZdcWestE->Fill(ZdcW_ene_horizontal+ZdcW_ene_vertical);

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

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

      TVector3 mom = femtoTrack->gMom();
      TVector3 dca = femtoTrack->gDCA(pVtx);

      hPt->Fill(mom.Pt());
      hPhi->Fill(mom.Phi());
      hEta->Fill(mom.Eta());

      hChi2->Fill(femtoTrack->chi2());
      hNhits->Fill(femtoTrack->nHits());
      hNhitsFit->Fill(femtoTrack->nHitsFit());
      hNhitsPoss->Fill(femtoTrack->nHitsPoss());

      hDedx->Fill(femtoTrack->dEdx());
      hTofM2->Fill(femtoTrack->massSqr());
      hCharge->Fill(femtoTrack->charge());
      
      hNsigEl->Fill(femtoTrack->nSigmaElectron());
      hNsigPi->Fill(femtoTrack->nSigmaPion());
      hNsigKa->Fill(femtoTrack->nSigmaKaon());
      hNsigPr->Fill(femtoTrack->nSigmaProton());

      hDCA->Fill(dca.Mag());

      // Simple single-track cut
      if( femtoTrack->gMom().Mag() < 0.1 || femtoTrack->gDCA(pVtx).Mag() > 3. ) {
	       continue;
      }

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

  fo->cd();
  
  // Vertex position
  hVtxX->Write();
  hVtxY->Write();
  hVtxZ->Write();
  hVpdZ->Write();
  // Energy from BBCs
  hBbcEastE->Write();
  hBbcWestE->Write();
  // Energy from ZDC-SMDs
  hZdcEastE->Write();
  hZdcWestE->Write();
  // Centrality (16 bins in total covering 0-80%)
  hCent16->Write();
  // refMults
  hRefMult->Write();
  hRefMult2->Write();
  hRefMultCorr->Write();
  hRefMultCorrW->Write();
  hNTofHit->Write();
  // Track information
  hPt->Write();
  hPhi->Write();
  hEta->Write();
  hChi2->Write();
  hNhits->Write();
  hNhitsFit->Write();
  hNhitsPoss->Write();
  hDedx->Write();
  hTofM2->Write();
  hNsigEl->Write();
  hNsigPi->Write();
  hNsigKa->Write();
  hNsigPr->Write();
  hDCA->Write();
  hCharge->Write();

  fo->Close();

  std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
	          << std::endl;
}