123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289 |
- /**
- * \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 "/mnt/pool/rhic/1/nigmatkulov/soft/StFemtoEvent/StFemtoDstReader.h"
- #include "/mnt/pool/rhic/1/nigmatkulov/soft/StFemtoEvent/StFemtoDst.h"
- #include "/mnt/pool/rhic/1/nigmatkulov/soft/StFemtoEvent/StFemtoEvent.h"
- #include "/mnt/pool/rhic/1/nigmatkulov/soft/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/1/nigmatkulov/soft/StFemtoEvent/libStFemtoDst)
- #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
- // 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);
- Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);
- // Used constants (global)
- const Int_t energy = 200.;
- // Used constants for the event cut
- const std::map<Float_t, Double_t> cutVtxZEnergy = {{7.7, 70.}, {11.5, 50.}, {19.6, 70}, {27., 70.}, {39., 40.}, {62., 40.}, {200., 30.}};
- const Double_t cutVtxR = 2.;
- const Double_t cutVpdVz = 3.;
- // Used constants for the track cut
- const std::map<Float_t, Double_t> cutDCA = {{7.7, 1.}, {11.5, 1.}, {19.6, 1.}, {27., 1.}, {39., 1.}, {62., 1.}, {200., 3.}};
- const Double_t cutEta = 1.;
- const Int_t cutNhits = 15;
- const Double_t cutNhitsPoss = 0.;
- const Double_t cutNhitsRatio = 0.51;
- const std::map<Float_t, Double_t> cutPtMin = {{7.7, 0.2}, {11.5, 0.2}, {19.6, 0.2}, {27., 0.2}, {39., 0.2}, {62., 0.2}, {200., 0.15}};
- const Double_t cutPtMax = 2.;
- const Double_t cutPMax = 10.;
- //_________________
- void FemtoDstAnalyzer(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",const Char_t *oFileName = "oTest.root") {
- std::cout << "Hi! Lets do some physics, Master!" << std::endl;
- #if ROOT_VERSION_CODE < ROOT_VERSION(6,0,0)
- gSystem->Load("/mnt/pool/rhic/1/nigmatkulov/soft/StFemtoEvent/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;
- // 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;
- }
- // Event selection
- if ( !isGoodEvent(event) ) continue;
- // Return primary vertex position
- TVector3 pVtx = event->primaryVertex();
- // Going through the inner BBC rings
- Double_t phiBBC_East, phiBBC_West;
- Double_t adcBBC_East, adcBBC_West;
- for (int iTile=0; iTile<16; iTile++)
- {
- // Get BBC East phi angle
- phiBBC_East = BBC_GetPhi(0,iTile);
- // Get BBC West phi angle
- phiBBC_West = BBC_GetPhi(1,iTile);
- // Get energy deposition from BBC East iTile module
- adcBBC_East = event->bbcAdcEast(iTile);
- // Get energy deposition from BBC West iTile module
- adcBBC_West = event->bbcAdcWest(iTile);
- }
- // 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 for EP determination in TPC
- if (!isGoodTrack(femtoTrack, energy, pVtx)) continue;
- //Track selection for flow calculation in TPC
- if (!isGoodTrackFlow(femtoTrack, energy, pVtx)) 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();
- std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
- << std::endl;
- }
- 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;
- }
- //--------------------------------------------------------------------------------------------------------------------------//
- //this function simply connects the gain values read in to the BBC azimuthal distribution
- //since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
- //I prefer assigning the angle between the tiles thus "greying out" the adcs.
- //Others have assigned all of the adc to one (exclusive) or the the other.
- Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId)
- {
- //float GetPhiInBBC(int eastWest, int bbcN) { //tileId=0 to 23
- const float Pi = TMath::Pi();
- const float Phi_div = Pi / 6;
- float bbc_phi = Phi_div;
- switch(tileId) {
- case 0: bbc_phi = 3.*Phi_div;
- break;
- case 1: bbc_phi = Phi_div;
- break;
- case 2: bbc_phi = -1.*Phi_div;
- break;
- case 3: bbc_phi = -3.*Phi_div;
- break;
- case 4: bbc_phi = -5.*Phi_div;
- break;
- case 5: bbc_phi = 5.*Phi_div;
- break;
- //case 6: bbc_phi= (mRndm.Rndm() > 0.5) ? 2.*Phi_div:4.*Phi_div; //tiles 7 and 9 are gained together we randomly assign the gain to one XOR the other
- case 6: bbc_phi = 3.*Phi_div;
- break;
- case 7: bbc_phi = 3.*Phi_div;
- break;
- case 8: bbc_phi = Phi_div;
- break;
- case 9: bbc_phi = 0.;
- break;
- case 10: bbc_phi = -1.*Phi_div;
- break;
- //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div; //tiles 13 and 15 are gained together
- case 11: bbc_phi = -3.*Phi_div;
- break;
- case 12: bbc_phi = -3.*Phi_div;
- break;
- case 13: bbc_phi = -5.*Phi_div;
- break;
- case 14: bbc_phi = Pi;
- break;
- case 15: bbc_phi = 5.*Phi_div;
- break;
- }
- //if we're looking at the east BBC we need to flip around x in the STAR coordinates,
- //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the
- if(0 == eastWest){
- if (bbc_phi > -0.001){ //this is not a >= since we are talking about finite adcs -- not to important
- bbc_phi = Pi - bbc_phi;
- }
- else {
- bbc_phi= -Pi - bbc_phi;
- }
- }
- if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
- if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;
- return bbc_phi;
- }
|