123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691 |
- /**
- * \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_EPQA(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
- const Char_t *outFile = "./oResTest.root",
- const Char_t *reCentFile = "",
- const Char_t *shiftFile = "")
- {
- 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;
- //Read ReCentering <Q> vectors
- TFile *fiReCent = new TFile(reCentFile, "read");
- TProfile2D *p_Qx_FULL_EP[fNharmonics][fArraySize];
- TProfile2D *p_Qy_FULL_EP[fNharmonics][fArraySize];
- TProfile2D *p_Qx_FULL_SP[fNharmonics][fArraySize];
- TProfile2D *p_Qy_FULL_SP[fNharmonics][fArraySize];
- TProfile2D *p_Qx_EAST_EP[fNharmonics][fArraySize][NEtaGaps];
- TProfile2D *p_Qy_EAST_EP[fNharmonics][fArraySize][NEtaGaps];
- TProfile2D *p_Qx_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
- TProfile2D *p_Qy_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
- TProfile2D *p_Qx_WEST_EP[fNharmonics][fArraySize][NEtaGaps];
- TProfile2D *p_Qy_WEST_EP[fNharmonics][fArraySize][NEtaGaps];
- TProfile2D *p_Qx_WEST_SP[fNharmonics][fArraySize][NEtaGaps];
- TProfile2D *p_Qy_WEST_SP[fNharmonics][fArraySize][NEtaGaps];
- for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
- {
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- p_Qx_FULL_EP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_EP%i", iHarm, iVtx));
- p_Qy_FULL_EP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_EP%i", iHarm, iVtx));
- p_Qx_FULL_SP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_SP%i", iHarm, iVtx));
- p_Qy_FULL_SP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_SP%i", iHarm, iVtx));
- for (int iEtaGap = 0; iEtaGap < NEtaGaps; iEtaGap++)
- {
- p_Qx_EAST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- p_Qy_EAST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- p_Qx_EAST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- p_Qy_EAST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- p_Qx_WEST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- p_Qy_WEST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- p_Qx_WEST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- p_Qy_WEST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- }
- }
- }
- //Read Shift coefficients
- TFile *fiShift = new TFile(shiftFile, "read");
- TProfile2D *p_Cos2_FULL_EP[fArraySize][NShiftOrderMax];
- TProfile2D *p_Sin2_FULL_EP[fArraySize][NShiftOrderMax];
- // TProfile2D *p_Cos2_FULL_SP[fArraySize][NShiftOrderMax];
- // TProfile2D *p_Sin2_FULL_SP[fArraySize][NShiftOrderMax];
- TProfile2D *p_Cos2_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin2_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Cos2_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin2_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Cos2_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin2_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Cos2_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin2_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Cos3_FULL_EP[fArraySize][NShiftOrderMax];
- TProfile2D *p_Sin3_FULL_EP[fArraySize][NShiftOrderMax];
- // TProfile2D *p_Cos3_FULL_SP[fArraySize][NShiftOrderMax];
- // TProfile2D *p_Sin3_FULL_SP[fArraySize][NShiftOrderMax];
- TProfile2D *p_Cos3_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin3_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Cos3_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin3_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Cos3_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin3_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Cos3_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- TProfile2D *p_Sin3_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
- {
- p_Cos2_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
- p_Sin2_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
- // p_Cos2_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
- // p_Sin2_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
- p_Cos3_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
- p_Sin3_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
- // p_Cos3_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
- // p_Sin3_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
- for (int iGap = 0; iGap < NEtaGaps; iGap++)
- {
- p_Cos2_EAST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Sin2_EAST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Cos2_EAST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Sin2_EAST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Cos2_WEST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Sin2_WEST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Cos2_WEST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Sin2_WEST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- p_Cos3_EAST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- p_Sin3_EAST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- p_Cos3_EAST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- p_Sin3_EAST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- p_Cos3_WEST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- p_Sin3_WEST_EP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- p_Cos3_WEST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- p_Sin3_WEST_SP[iVtx][iShift][iGap] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- }
- }
- }
- //Manage output TProfiles
- TH2D *hPsi2EPRaw[NEtaGaps];
- TH2D *hPsi2EPRec[NEtaGaps];
- TH2D *hPsi2EPFlt[NEtaGaps];
- TH2D *hPsi3EPRaw[NEtaGaps];
- TH2D *hPsi3EPRec[NEtaGaps];
- TH2D *hPsi3EPFlt[NEtaGaps];
-
- for (int iGap=0;iGap<NEtaGaps;iGap++)
- {
- hPsi2EPRaw[iGap] = new TH2D(Form("hPsi2EPRaw%i",iGap),Form("hPsi2EPRaw%i",iGap),500,-5./2.,5./2.,9,-0.5,8.5);
- hPsi2EPRec[iGap] = new TH2D(Form("hPsi2EPRec%i",iGap),Form("hPsi2EPRec%i",iGap),500,-5./2.,5./2.,9,-0.5,8.5);
- hPsi2EPFlt[iGap] = new TH2D(Form("hPsi2EPFlt%i",iGap),Form("hPsi2EPFlt%i",iGap),500,-5./2.,5./2.,9,-0.5,8.5);
- hPsi3EPRaw[iGap] = new TH2D(Form("hPsi3EPRaw%i",iGap),Form("hPsi2EPRaw%i",iGap),500,-5./3.,5./3.,9,-0.5,8.5);
- hPsi3EPRec[iGap] = new TH2D(Form("hPsi3EPRec%i",iGap),Form("hPsi2EPRec%i",iGap),500,-5./3.,5./3.,9,-0.5,8.5);
- hPsi3EPFlt[iGap] = new TH2D(Form("hPsi3EPFlt%i",iGap),Form("hPsi2EPFlt%i",iGap),500,-5./3.,5./3.,9,-0.5,8.5);
- }
- //Counters
- Int_t fUsedTracks = 0;
- Int_t fQcountFull = 0;
- Int_t fQcountFullEast = 0;
- Int_t fQcountFullWest = 0;
- Int_t fQcountWest[NEtaGaps];
- Int_t fQcountEast[NEtaGaps];
- //Q-vectors
- TVector2 fQv2FullEP;
- TVector2 fQv2FullSP;
- TVector2 fQv3FullEP;
- TVector2 fQv3FullSP;
- TVector2 fQv2FullEPRaw;
- TVector2 fQv3FullEPRaw;
- TVector2 fQv2WestEP[NEtaGaps];
- TVector2 fQv2EastEP[NEtaGaps];
- TVector2 fQv2WestSP[NEtaGaps];
- TVector2 fQv2EastSP[NEtaGaps];
- TVector2 fQv3WestEP[NEtaGaps];
- TVector2 fQv3EastEP[NEtaGaps];
- TVector2 fQv3WestSP[NEtaGaps];
- TVector2 fQv3EastSP[NEtaGaps];
- TVector2 QvTMP;
- Double_t weight;
- Int_t VtxSign;
- Double_t psiShiftedTMP;
- Double_t PsiSin;
- Double_t PsiCos;
- Long64_t goodEventCounter = 0;
- Double_t Psi2EastEP[NEtaGaps], Psi3EastEP[NEtaGaps];
- Double_t Psi2WestEP[NEtaGaps], Psi3WestEP[NEtaGaps];
- Double_t Psi_ReCenter, dPsi;
- Int_t bin;
- Double_t res2EP[NEtaGaps], res3EP[NEtaGaps];
- // 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();
- // Track analysis
- Int_t nTracks = dst->numberOfTracks();
-
- //Counters
- fUsedTracks = 0;
- fQcountFull = 0;
- fQcountFullEast = 0;
- fQcountFullWest = 0;
- psiShiftedTMP = 0.;
- PsiSin = 0.;
- PsiCos = 0.;
- fQv2FullEP.Set(0.,0.);
- fQv2FullEPRaw.Set(0.,0.);
- fQv2FullSP.Set(0.,0.);
- fQv3FullEP.Set(0.,0.);
- fQv3FullEPRaw.Set(0.,0.);
- fQv3FullSP.Set(0.,0.);
- QvTMP.Set(0.,0.);
- for (int i=0; i<NEtaGaps; i++)
- {
- fQcountWest[i] = 0;
- fQcountEast[i] = 0;
- fQv2WestEP[i].Set(0.,0.);
- fQv2EastEP[i].Set(0.,0.);
- fQv2WestSP[i].Set(0.,0.);
- fQv2EastSP[i].Set(0.,0.);
- fQv3WestEP[i].Set(0.,0.);
- fQv3EastEP[i].Set(0.,0.);
- fQv3WestSP[i].Set(0.,0.);
- fQv3EastSP[i].Set(0.,0.);
- }
- // Vertex sign
- VtxSign = (event->primaryVertex().Z() > 0) ? 0 : 1;
- // 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 (!isGoodTrack(femtoTrack, energy, pVtx)) continue;
- fUsedTracks++;
-
- //Fill recentered Q-vectors for all TPC (FULL)
- weight = GetWeight(femtoTrack);
-
- QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign]->GetBinContent(p_Qx_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_FULL_EP[0][VtxSign]->GetBinContent(p_Qy_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv2FullEPRaw += weight*(CalcQvector(femtoTrack,2));
- fQv2FullEP += weight*(CalcQvector(femtoTrack,2) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_FULL_EP[1][VtxSign]->GetBinContent(p_Qx_FULL_EP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_FULL_EP[1][VtxSign]->GetBinContent(p_Qy_FULL_EP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv3FullEPRaw += weight*(CalcQvector(femtoTrack,3));
- fQv3FullEP += weight*(CalcQvector(femtoTrack,3) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_FULL_SP[0][VtxSign]->GetBinContent(p_Qx_FULL_SP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_FULL_SP[0][VtxSign]->GetBinContent(p_Qy_FULL_SP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
-
- fQv2FullSP += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_FULL_SP[1][VtxSign]->GetBinContent(p_Qx_FULL_SP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_FULL_SP[1][VtxSign]->GetBinContent(p_Qy_FULL_SP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv3FullSP += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
- QvTMP.Set(0.,0.);
-
- fQcountFull++;
- if (femtoTrack->eta() >= 0.)
- {
- fQcountFullWest++;
- }
- else
- {
- fQcountFullEast++;
- }
- for (int iGap=0; iGap<NEtaGaps; iGap++)
- {
- //fill Qvectors from EAST arm
- if (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < cutEtaGap[iGap])
- {
- QvTMP.SetX(p_Qx_EAST_EP[0][VtxSign][iGap]->GetBinContent(p_Qx_EAST_EP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_EAST_EP[0][VtxSign][iGap]->GetBinContent(p_Qy_EAST_EP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv2EastEP[iGap] += weight*(CalcQvector(femtoTrack,2) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_EAST_EP[1][VtxSign][iGap]->GetBinContent(p_Qx_EAST_EP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_EAST_EP[1][VtxSign][iGap]->GetBinContent(p_Qy_EAST_EP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv3EastEP[iGap] += weight*(CalcQvector(femtoTrack,3) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_EAST_SP[0][VtxSign][iGap]->GetBinContent(p_Qx_EAST_SP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_EAST_SP[0][VtxSign][iGap]->GetBinContent(p_Qy_EAST_SP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
-
- fQv2EastSP[iGap] += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_EAST_SP[1][VtxSign][iGap]->GetBinContent(p_Qx_EAST_SP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_EAST_SP[1][VtxSign][iGap]->GetBinContent(p_Qy_EAST_SP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv2EastSP[iGap] += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
- QvTMP.Set(0.,0.);
- fQcountEast[iGap]++;
- }
- //fill Qvectors from WEST arm
- if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > cutEtaGap[iGap])
- {
- QvTMP.SetX(p_Qx_WEST_EP[0][VtxSign][iGap]->GetBinContent(p_Qx_WEST_EP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_WEST_EP[0][VtxSign][iGap]->GetBinContent(p_Qy_WEST_EP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv2WestEP[iGap] += weight*(CalcQvector(femtoTrack,2) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_WEST_EP[1][VtxSign][iGap]->GetBinContent(p_Qx_WEST_EP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_WEST_EP[1][VtxSign][iGap]->GetBinContent(p_Qy_WEST_EP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv3WestEP[iGap] += weight*(CalcQvector(femtoTrack,3) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_WEST_SP[0][VtxSign][iGap]->GetBinContent(p_Qx_WEST_SP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_WEST_SP[0][VtxSign][iGap]->GetBinContent(p_Qy_WEST_SP[0][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
-
- fQv2WestSP[iGap] += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
- QvTMP.Set(0.,0.);
- QvTMP.SetX(p_Qx_WEST_SP[1][VtxSign][iGap]->GetBinContent(p_Qx_WEST_SP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- QvTMP.SetY(p_Qy_WEST_SP[1][VtxSign][iGap]->GetBinContent(p_Qy_WEST_SP[1][VtxSign][iGap]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
- fQv2WestSP[iGap] += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
- QvTMP.Set(0.,0.);
- fQcountWest[iGap]++;
- }
- }
- // Check if track has TOF signal
- if (femtoTrack->isTofTrack())
- {
- } //if( isTofTrack() )
- } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
- hPsi2EPRaw[0]->Fill((Double_t)TMath::ATan2(fQv2FullEPRaw.Y(),fQv2FullEPRaw.X())/2.,(Double_t)event->cent9());
- hPsi2EPRec[0]->Fill((Double_t)TMath::ATan2(fQv2FullEP.Y(),fQv2FullEP.X())/2.,(Double_t)event->cent9());
- hPsi3EPRaw[0]->Fill((Double_t)TMath::ATan2(fQv3FullEPRaw.Y(),fQv3FullEPRaw.X())/3.,(Double_t)event->cent9());
- hPsi3EPRec[0]->Fill((Double_t)TMath::ATan2(fQv3FullEP.Y(),fQv3FullEP.X())/3.,(Double_t)event->cent9());
- //Getting PsiEP with recentering + shift correction
- Psi_ReCenter=0.;
- Double_t mCos, mSin;
- Double_t Psi2EPFlat, Psi3EPFlat, Psi2EP, Psi3EP;
- Psi2EP = TMath::ATan2(fQv2FullEP.Y(),fQv2FullEP.X())/2.;
- Psi3EP = TMath::ATan2(fQv3FullEP.Y(),fQv3FullEP.X())/3.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mCos = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);
- bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mSin = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);
- dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi2EP));
- }
- Psi2EPFlat = AngleShift(Psi2EP+dPsi,2.);
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- bin = p_Cos3_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mCos = p_Cos3_FULL_EP[VtxSign][iK]->GetBinContent(bin);
- bin = p_Cos3_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mSin = p_Cos3_FULL_EP[VtxSign][iK]->GetBinContent(bin);
- dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff3.at(iK)*Psi3EP)+
- mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi3EP));
- }
- Psi3EPFlat = AngleShift(Psi3EP+dPsi,3.);
- hPsi2EPFlt[0]->Fill(Psi2EPFlat,(Double_t)event->cent9());
- hPsi3EPFlt[0]->Fill(Psi3EPFlat,(Double_t)event->cent9());
- for (int iGap=0; iGap<NEtaGaps; iGap++)
- {
- if (fQcountEast[iGap] > cutNcountQvGap &&
- fQcountWest[iGap] > cutNcountQvGap)
- {
- Psi2EastEP[iGap]=0.;
- Psi3EastEP[iGap]=0.;
- Psi2WestEP[iGap]=0.;
- Psi3WestEP[iGap]=0.;
- //EAST 2-nd order
- Psi_ReCenter = TMath::ATan2(fQv2EastEP[iGap].Y(),fQv2EastEP[iGap].X())/2.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- bin = p_Cos2_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mCos = p_Cos2_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- bin = p_Cos2_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mSin = p_Cos2_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
- }
- Psi2EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
- //WEST 2-nd order
- Psi_ReCenter = TMath::ATan2(fQv2WestEP[iGap].Y(),fQv2WestEP[iGap].X())/2.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- bin = p_Cos2_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mCos = p_Cos2_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- bin = p_Cos2_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mSin = p_Cos2_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
- }
- Psi2WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
- //EAST 3-nd order
- Psi_ReCenter = TMath::ATan2(fQv3EastEP[iGap].Y(),fQv3EastEP[iGap].X())/3.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- bin = p_Cos3_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mCos = p_Cos3_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- bin = p_Cos3_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mSin = p_Cos3_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
- }
- Psi3EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
- //WEST 3-nd order
- Psi_ReCenter = TMath::ATan2(fQv3WestEP[iGap].Y(),fQv3WestEP[iGap].X())/3.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- bin = p_Cos3_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mCos = p_Cos3_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- bin = p_Cos3_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
- mSin = p_Cos3_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
- dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
- }
- Psi3WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
- }
- }
- } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
- femtoReader->Finish();
- TFile *output = new TFile(outFile,"recreate");
- for (int iGap=0;iGap<NEtaGaps;iGap++)
- {
- hPsi2EPRaw[iGap]->Write();
- hPsi2EPRec[iGap]->Write();
- hPsi2EPFlt[iGap]->Write();
- hPsi3EPRaw[iGap]->Write();
- hPsi3EPRec[iGap]->Write();
- hPsi3EPFlt[iGap]->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;
- }
- 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;
- }
|