StFlowMephi/FemtoDstAnalyzer_BBC_Gain.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"
#include "TRandom3.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);
Int_t GetVzBin(Double_t vtxZ);

//_________________
void FemtoDstAnalyzer_BBC_Gain(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 *BBCgain_East[fArraySize];
  TProfile2D *BBCgain_West[fArraySize];

  for (int iVtx=0; iVtx<fArraySize; iVtx++)
  {
    BBCgain_East[iVtx] = new TProfile2D(Form("p_gainBBC_East_Vtx%i",iVtx),Form("p_gainBBC_East_Vtx%i",iVtx),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
    BBCgain_West[iVtx] = new TProfile2D(Form("p_gainBBC_West_Vtx%i",iVtx),Form("p_gainBBC_West_Vtx%i",iVtx),fNBins,runId_min-0.5,runId_max+0.5,16,-0.5,15.5);
  }


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

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

    for (int iModule=0; iModule<24; iModule++)
    {
      BBCgain_East[VtxSign]->Fill(event->runId(),iModule,event->bbcAdcEast(iModule));
      BBCgain_West[VtxSign]->Fill(event->runId(),iModule,event->bbcAdcWest(iModule));
    }

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

  femtoReader->Finish();

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

  for (int iVtx=0; iVtx<fArraySize; iVtx++)
  {
    BBCgain_East[iVtx]->Write();
    BBCgain_West[iVtx]->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;
}

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