|
@@ -0,0 +1,287 @@
|
|
|
+/**
|
|
|
+ * \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 "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDstReader.h"
|
|
|
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDst.h"
|
|
|
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoEvent.h"
|
|
|
+#include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/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);
|
|
|
+Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
|
|
|
+Double_t GetWeight(StFemtoTrack *const &track);
|
|
|
+TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
|
|
|
+Double_t AngleShift(Double_t Psi, Double_t order);
|
|
|
+
|
|
|
+//_________________
|
|
|
+void FemtoDstAnalyzer_PID(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
|
|
|
+ const Char_t *outFile = "./oPIDTest.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;
|
|
|
+
|
|
|
+
|
|
|
+ const int NptBins = 24;
|
|
|
+ const double ptBin [NptBins+1] = {0.,0.2,0.4,0.6,0.8,1.,1.2,1.4,1.6,1.8,2.,2.2,2.4,2.6,2.8,3.,3.2,3.4,3.6,3.8,4.,4.5,5.,5.5,6.};
|
|
|
+ TH2F *hNsigPionMSqr[NptBins];
|
|
|
+ TH2F *hNsigKaonMSqr[NptBins];
|
|
|
+ TH2F *hNsigProtonMSqr[NptBins];
|
|
|
+
|
|
|
+ //Initialization
|
|
|
+ for (int i=0; i<NptBins; i++)
|
|
|
+ {
|
|
|
+ hNsigPionMSqr[i] = new TH2F(Form("hNsigPionMSqr%i",i) ,Form("n#sigma_{#pi} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{#pi};m^{2} [GeV/c^{2}]^{2}",ptBin[i]+cutPtMin.at(energy),ptBin[i+1]+cutPtMin.at(energy)),1400,-35.,35.,1400,-2.,5.);
|
|
|
+ hNsigKaonMSqr[i] = new TH2F(Form("hNsigKaonMSqr%i",i) ,Form("n#sigma_{K} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{K};m^{2} [GeV/c^{2}]^{2}",ptBin[i]+cutPtMin.at(energy),ptBin[i+1]+cutPtMin.at(energy)),1400,-35.,35.,1400,-2.,5.);
|
|
|
+ hNsigProtonMSqr[i] = new TH2F(Form("hNsigProtonMSqr%i",i) ,Form("n#sigma_{p} vs M^{2} for %.2f < p_{T} < %.2f GeV/c;n#sigma_{p};m^{2} [GeV/c^{2}]^{2}",ptBin[i]+cutPtMin.at(energy),ptBin[i+1]+cutPtMin.at(energy)),1400,-35.,35.,1400,-2.,5.);
|
|
|
+ }
|
|
|
+
|
|
|
+ // 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;
|
|
|
+ }
|
|
|
+
|
|
|
+ TVector3 pVtx = event->primaryVertex();
|
|
|
+
|
|
|
+ // Event selection
|
|
|
+ if ( !isGoodEvent(event) ) continue;
|
|
|
+
|
|
|
+ // 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;
|
|
|
+
|
|
|
+ //Track selection
|
|
|
+ if (!isGoodTrackFlow(femtoTrack, energy, pVtx)) continue;
|
|
|
+
|
|
|
+ // Determine pt bin for PID plots
|
|
|
+ int i_pt = -1;
|
|
|
+ for (int i=0;i<NptBins;i++)
|
|
|
+ {
|
|
|
+ if (femtoTrack->pt() >= (ptBin[i] + cutPtMin.at(energy)) && femtoTrack->pt() < (ptBin[i+1] + cutPtMin.at(energy)))
|
|
|
+ {
|
|
|
+ i_pt = i;
|
|
|
+ }
|
|
|
+ }
|
|
|
+ if (i_pt == -1) continue;
|
|
|
+
|
|
|
+ // Check if track has TOF signal
|
|
|
+ if (femtoTrack->isTofTrack())
|
|
|
+ {
|
|
|
+ if (femtoTrack->gDCA(pVtx).Mag() >= 1.) continue;
|
|
|
+
|
|
|
+ hNsigPionMSqr[i_pt]->Fill(femtoTrack->nSigmaPion(),femtoTrack->massSqr());
|
|
|
+ hNsigKaonMSqr[i_pt]->Fill(femtoTrack->nSigmaKaon(),femtoTrack->massSqr());
|
|
|
+ hNsigProtonMSqr[i_pt]->Fill(femtoTrack->nSigmaProton(),femtoTrack->massSqr());
|
|
|
+
|
|
|
+ } //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 i=0; i<NptBins; i++)
|
|
|
+ {
|
|
|
+ hNsigPionMSqr[i] -> Write();
|
|
|
+ hNsigKaonMSqr[i] -> Write();
|
|
|
+ hNsigProtonMSqr[i] -> 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;
|
|
|
+}
|
|
|
+
|
|
|
+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;
|
|
|
+}
|
|
|
+
|
|
|
+Bool_t isGoodTrackFlow(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;
|
|
|
+}
|
|
|
+
|
|
|
+Double_t AngleShift(Double_t Psi, Double_t order)
|
|
|
+{
|
|
|
+ Double_t PsiCorr = Psi;
|
|
|
+ if (PsiCorr > TMath::Pi()/order)
|
|
|
+ {
|
|
|
+ PsiCorr = Psi - 2.*TMath::Pi()/order;
|
|
|
+ }
|
|
|
+ if (PsiCorr < -1.*TMath::Pi()/order)
|
|
|
+ {
|
|
|
+ PsiCorr = Psi + 2.*TMath::Pi()/order;
|
|
|
+ }
|
|
|
+ return PsiCorr;
|
|
|
+}
|