123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375 |
- /**
- * \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);
- Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip);
- //_________________
- void FemtoDstAnalyzer_ZDC_ShiftCorr(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
- const Char_t *outFile = "./oReCenteringTest.root",
- const Char_t *recZDC = "")
- {
- 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
- //Manage recentering
- TProfile2D *p_ZDC_Qx_Full_EP[fNharmonics+1][fArraySize];
- TProfile2D *p_ZDC_Qy_Full_EP[fNharmonics+1][fArraySize];
- TProfile2D *p_ZDC_Qx_East_EP[fNharmonics+1][fArraySize];
- TProfile2D *p_ZDC_Qy_East_EP[fNharmonics+1][fArraySize];
- TProfile2D *p_ZDC_Qx_West_EP[fNharmonics+1][fArraySize];
- TProfile2D *p_ZDC_Qy_West_EP[fNharmonics+1][fArraySize];
- TFile *fiReCentZDC = new TFile(recZDC, "read");
- fiReCentZDC->cd();
- for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
- {
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- p_ZDC_Qx_Full_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
- p_ZDC_Qy_Full_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
- p_ZDC_Qx_East_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
- p_ZDC_Qy_East_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
- p_ZDC_Qx_West_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
- p_ZDC_Qy_West_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
- }
- }
- //Manage output
- TProfile2D *p_ZDC_Cos_Full_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
- TProfile2D *p_ZDC_Sin_Full_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
- TProfile2D *p_ZDC_Cos_East_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
- TProfile2D *p_ZDC_Sin_East_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
- TProfile2D *p_ZDC_Cos_West_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
- TProfile2D *p_ZDC_Sin_West_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
-
- for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
- {
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
- {
- p_ZDC_Cos_Full_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- p_ZDC_Sin_Full_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- p_ZDC_Cos_East_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- p_ZDC_Sin_East_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- p_ZDC_Cos_West_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- p_ZDC_Sin_West_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- }
- }
- }
- 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;
- Int_t VtxSign;
- // 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;
- VtxSign = GetVzBin(event->primaryVertex().Z());
- if (VtxSign == -1)
- continue;
- if (event->zdcSumAdcEast() < 1.) continue;
- if (event->zdcSumAdcWest() < 1.) continue;
- Double_t qxEast, qyEast, adcxEast, adcyEast;
- Double_t qxWest, qyWest, adcxWest, adcyWest;
- Double_t qxFull, qyFull, adcxFull, adcyFull;
- qxEast = 0.;
- qyEast = 0.;
- adcxEast = 0.;
- adcyEast = 0.;
- qxWest = 0.;
- qyWest = 0.;
- adcxWest = 0.;
- adcyWest = 0.;
- qxFull = 0.;
- qyFull = 0.;
- adcxFull = 0.;
- adcyFull = 0.;
- Double_t qxRecEast, qyRecEast;
- Double_t qxRecWest, qyRecWest;
- Double_t qxRecFull, qyRecFull;
- // Vertical ZDC-SMD
- for (int iStrip=0; iStrip<fNZdcSmdStripsVertical; iStrip++)
- {
- qxEast += event->zdcSmdEastVertical(iStrip) * GetZDCPosition(0,0,iStrip);
- adcxEast += event->zdcSmdEastVertical(iStrip);
- qxWest += event->zdcSmdWestVertical(iStrip) * GetZDCPosition(1,0,iStrip);
- adcxWest += event->zdcSmdWestVertical(iStrip);
- }
- // Horizontal ZDC-SMD
- for (int iStrip=0; iStrip<fNZdcSmdStripsHorizontal; iStrip++)
- {
- qyEast += event->zdcSmdEastHorizontal(iStrip) * GetZDCPosition(0,1,iStrip);
- adcyEast += event->zdcSmdEastHorizontal(iStrip);
- qyWest += event->zdcSmdWestHorizontal(iStrip) * GetZDCPosition(1,1,iStrip);
- adcyWest += event->zdcSmdWestHorizontal(iStrip);
- }
- if (adcxEast > 0) qxEast /= adcxEast;
- else qxEast = -999.;
- if (adcyEast > 0) qyEast /= adcyEast;
- else qyEast = -999.;
- if (adcxWest > 0) qxWest /= adcxWest;
- else qxWest = -999.;
- if (adcyWest > 0) qyWest /= adcyWest;
- else qyWest = -999.;
- // Full ZDC-SMD Q-vector
- if (qxEast == -999. || qxWest == -999.) qxFull = -999.;
- else qxFull = qxEast - qxWest;
- if (qyEast == -999. || qyWest == -999.) qyFull = -999.;
- else qyFull = qyEast - qyWest;
- if (qxEast != -999.) qxRecEast = qxEast - (Double_t)p_ZDC_Qx_East_EP[0][VtxSign]->GetBinContent(p_ZDC_Qx_East_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
- if (qyEast != -999.) qyRecEast = qyEast - (Double_t)p_ZDC_Qy_East_EP[0][VtxSign]->GetBinContent(p_ZDC_Qy_East_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
- if (qxWest != -999.) qxRecWest = qxWest - (Double_t)p_ZDC_Qx_West_EP[0][VtxSign]->GetBinContent(p_ZDC_Qx_West_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
- if (qyWest != -999.) qyRecWest = qyWest - (Double_t)p_ZDC_Qy_West_EP[0][VtxSign]->GetBinContent(p_ZDC_Qy_West_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
- if (qxFull != -999.) qxRecFull = qxFull - (Double_t)p_ZDC_Qx_Full_EP[0][VtxSign]->GetBinContent(p_ZDC_Qx_Full_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
- if (qyFull != -999.) qyRecFull = qyFull - (Double_t)p_ZDC_Qy_Full_EP[0][VtxSign]->GetBinContent(p_ZDC_Qy_Full_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
- if (qxEast == -999.) qxRecEast = qxEast;
- if (qyEast == -999.) qyRecEast = qyEast;
- if (qxWest == -999.) qxRecWest = qxWest;
- if (qyWest == -999.) qyRecWest = qyWest;
- if (qxFull == -999.) qxRecFull = qxFull;
- if (qyFull == -999.) qyRecFull = qyFull;
- Double_t PsiSin = 0., PsiCos = 0., psiShiftedTMP = 0.;
- TVector2 qRecEast[fNharmonics+1], qRecWest[fNharmonics+1], qRecFull[fNharmonics+1];
- for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
- {
- qRecEast[iHarm].Set(qxRecEast,qyRecEast);
- qRecWest[iHarm].Set(qxRecWest,qyRecWest);
- qRecFull[iHarm].Set(qxRecFull,qyRecFull);
- for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
- {
- //Full ZDC
- if (qRecFull[iHarm].Mod())
- {
- psiShiftedTMP = qRecFull[iHarm].Phi() / (iHarm + 1.);
- if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
- }
- PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
- PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
-
- p_ZDC_Cos_Full_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
- p_ZDC_Sin_Full_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);
- //East ZDC
- if (qRecEast[iHarm].Mod())
- {
- psiShiftedTMP = qRecEast[iHarm].Phi() / (iHarm + 1.);
- if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
- }
- PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
- PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
-
- p_ZDC_Cos_East_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
- p_ZDC_Sin_East_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);
- //West ZDC
- if (qRecWest[iHarm].Mod())
- {
- psiShiftedTMP = qRecWest[iHarm].Phi() / (iHarm + 1.);
- if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
- }
- PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
- PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
-
- p_ZDC_Cos_West_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
- p_ZDC_Sin_West_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);
- }
- }
- } // end of the event loop
-
- femtoReader->Finish();
- TFile *output = new TFile(outFile, "recreate");
- for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
- {
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
- {
- p_ZDC_Cos_Full_EP[iHarm][iVtx][iShift]->Write();
- p_ZDC_Sin_Full_EP[iHarm][iVtx][iShift]->Write();
- p_ZDC_Cos_East_EP[iHarm][iVtx][iShift]->Write();
- p_ZDC_Sin_East_EP[iHarm][iVtx][iShift]->Write();
- p_ZDC_Cos_West_EP[iHarm][iVtx][iShift]->Write();
- p_ZDC_Sin_West_EP[iHarm][iVtx][iShift]->Write();
- }
- }
- }
- output->Close();
- std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
- << 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;
- }
- Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip)
- // Get position of each slat;strip starts from 0
- {
- Double_t zdcsmd_x[7] = {0.5,2,3.5,5,6.5,8,9.5};
- Double_t zdcsmd_y[8] = {1.25,3.25,5.25,7.25,9.25,11.25,13.25,15.25};
- Double_t mZDCSMDCenterex = 4.72466;
- Double_t mZDCSMDCenterey = 5.53629;
- Double_t mZDCSMDCenterwx = 4.39604;
- Double_t mZDCSMDCenterwy = 5.19968;
- if(eastwest==0 && verthori==0) return zdcsmd_x[strip]-mZDCSMDCenterex;
- if(eastwest==1 && verthori==0) return mZDCSMDCenterwx-zdcsmd_x[strip];
- if(eastwest==0 && verthori==1) return zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterey;
- if(eastwest==1 && verthori==1) return zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterwy;
- return -999.;
- }
|