FemtoDstAnalyzer_EPQA.C 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691
  1. /**
  2. * \brief Example of how to read a file (list of files) using StFemtoEvent classes
  3. *
  4. * RunFemtoDstAnalyzer.C is an example of reading FemtoDst format.
  5. * One can use either FemtoDst file or a list of femtoDst files (inFile.lis or
  6. * inFile.list) as an input, and preform physics analysis
  7. *
  8. * \author Grigory Nigmatkulov
  9. * \date May 29, 2018
  10. */
  11. // This is needed for calling standalone classes
  12. #define _VANILLA_ROOT_
  13. // C++ headers
  14. #include <iostream>
  15. #include <vector>
  16. #include <map>
  17. // ROOT headers
  18. #include "TROOT.h"
  19. #include "TFile.h"
  20. #include "TChain.h"
  21. #include "TVector2.h"
  22. #include "TVector3.h"
  23. #include "TTree.h"
  24. #include "TSystem.h"
  25. #include "TH1.h"
  26. #include "TH2.h"
  27. #include "TMath.h"
  28. #include "TProfile2D.h"
  29. #include "TStopwatch.h"
  30. // FemtoDst headers
  31. #include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDstReader.h"
  32. #include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoDst.h"
  33. #include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoEvent.h"
  34. #include "/mnt/pool/rhic/4/parfenovpeter/STAR/StFemtoEvent/StFemtoTrack.h"
  35. // Load libraries (for ROOT_VERSTION_CODE >= 393215)
  36. #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
  37. R__LOAD_LIBRARY(/mnt/pool/rhic/4/parfenovpeter/STAR/build/libStFemtoDst.so)
  38. #endif
  39. #include "Constants.h"
  40. // inFile - is a name of name.FemtoDst.root file or a name
  41. // of a name.lis(t) files, that contains a list of
  42. // name1.FemtoDst.root, name2.FemtoDst.root, ... files
  43. //Used functions (see them below)
  44. Bool_t isGoodEvent(StFemtoEvent *const &event);
  45. Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
  46. Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
  47. Double_t GetWeight(StFemtoTrack *const &track);
  48. TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
  49. Double_t AngleShift(Double_t Psi, Double_t order);
  50. //_________________
  51. void FemtoDstAnalyzer_EPQA(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
  52. const Char_t *outFile = "./oResTest.root",
  53. const Char_t *reCentFile = "",
  54. const Char_t *shiftFile = "")
  55. {
  56. std::cout << "Hi! Lets do some physics, Master!" << std::endl;
  57. TStopwatch timer;
  58. timer.Start();
  59. #if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
  60. gSystem->Load("../libStFemtoDst.so");
  61. #endif
  62. StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
  63. femtoReader->Init();
  64. // This is a way if you want to spead up IO
  65. std::cout << "Explicit read status for some branches" << std::endl;
  66. femtoReader->SetStatus("*", 0);
  67. femtoReader->SetStatus("Event", 1);
  68. femtoReader->SetStatus("Track", 1);
  69. std::cout << "Status has been set" << std::endl;
  70. std::cout << "Now I know what to read, Master!" << std::endl;
  71. if (!femtoReader->chain())
  72. {
  73. std::cout << "No chain has been found." << std::endl;
  74. }
  75. Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  76. std::cout << "eventsInTree: " << eventsInTree << std::endl;
  77. Long64_t events2read = femtoReader->chain()->GetEntries();
  78. std::cout << "Number of events to read: " << events2read << std::endl;
  79. //Read ReCentering <Q> vectors
  80. TFile *fiReCent = new TFile(reCentFile, "read");
  81. TProfile2D *p_Qx_FULL_EP[fNharmonics][fArraySize];
  82. TProfile2D *p_Qy_FULL_EP[fNharmonics][fArraySize];
  83. TProfile2D *p_Qx_FULL_SP[fNharmonics][fArraySize];
  84. TProfile2D *p_Qy_FULL_SP[fNharmonics][fArraySize];
  85. TProfile2D *p_Qx_EAST_EP[fNharmonics][fArraySize][NEtaGaps];
  86. TProfile2D *p_Qy_EAST_EP[fNharmonics][fArraySize][NEtaGaps];
  87. TProfile2D *p_Qx_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
  88. TProfile2D *p_Qy_EAST_SP[fNharmonics][fArraySize][NEtaGaps];
  89. TProfile2D *p_Qx_WEST_EP[fNharmonics][fArraySize][NEtaGaps];
  90. TProfile2D *p_Qy_WEST_EP[fNharmonics][fArraySize][NEtaGaps];
  91. TProfile2D *p_Qx_WEST_SP[fNharmonics][fArraySize][NEtaGaps];
  92. TProfile2D *p_Qy_WEST_SP[fNharmonics][fArraySize][NEtaGaps];
  93. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  94. {
  95. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  96. {
  97. p_Qx_FULL_EP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_EP%i", iHarm, iVtx));
  98. p_Qy_FULL_EP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_EP%i", iHarm, iVtx));
  99. p_Qx_FULL_SP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_SP%i", iHarm, iVtx));
  100. p_Qy_FULL_SP[iHarm][iVtx] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_SP%i", iHarm, iVtx));
  101. for (int iEtaGap = 0; iEtaGap < NEtaGaps; iEtaGap++)
  102. {
  103. p_Qx_EAST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  104. p_Qy_EAST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  105. p_Qx_EAST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  106. p_Qy_EAST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  107. p_Qx_WEST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  108. p_Qy_WEST_EP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  109. p_Qx_WEST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  110. p_Qy_WEST_SP[iHarm][iVtx][iEtaGap] = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  111. }
  112. }
  113. }
  114. //Read Shift coefficients
  115. TFile *fiShift = new TFile(shiftFile, "read");
  116. TProfile2D *p_Cos2_FULL_EP[fArraySize][NShiftOrderMax];
  117. TProfile2D *p_Sin2_FULL_EP[fArraySize][NShiftOrderMax];
  118. // TProfile2D *p_Cos2_FULL_SP[fArraySize][NShiftOrderMax];
  119. // TProfile2D *p_Sin2_FULL_SP[fArraySize][NShiftOrderMax];
  120. TProfile2D *p_Cos2_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  121. TProfile2D *p_Sin2_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  122. TProfile2D *p_Cos2_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  123. TProfile2D *p_Sin2_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  124. TProfile2D *p_Cos2_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  125. TProfile2D *p_Sin2_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  126. TProfile2D *p_Cos2_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  127. TProfile2D *p_Sin2_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  128. TProfile2D *p_Cos3_FULL_EP[fArraySize][NShiftOrderMax];
  129. TProfile2D *p_Sin3_FULL_EP[fArraySize][NShiftOrderMax];
  130. // TProfile2D *p_Cos3_FULL_SP[fArraySize][NShiftOrderMax];
  131. // TProfile2D *p_Sin3_FULL_SP[fArraySize][NShiftOrderMax];
  132. TProfile2D *p_Cos3_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  133. TProfile2D *p_Sin3_EAST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  134. TProfile2D *p_Cos3_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  135. TProfile2D *p_Sin3_EAST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  136. TProfile2D *p_Cos3_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  137. TProfile2D *p_Sin3_WEST_EP[fArraySize][NEtaGaps][NShiftOrderMax];
  138. TProfile2D *p_Cos3_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  139. TProfile2D *p_Sin3_WEST_SP[fArraySize][NEtaGaps][NShiftOrderMax];
  140. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  141. {
  142. for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
  143. {
  144. p_Cos2_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
  145. p_Sin2_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
  146. // p_Cos2_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
  147. // p_Sin2_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
  148. p_Cos3_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
  149. p_Sin3_FULL_EP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
  150. // p_Cos3_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
  151. // p_Sin3_FULL_SP[iVtx][iShift] = (TProfile2D*) fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
  152. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  153. {
  154. 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));
  155. 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));
  156. 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));
  157. 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));
  158. 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));
  159. 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));
  160. 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));
  161. 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));
  162. 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));
  163. 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));
  164. 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));
  165. 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));
  166. 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));
  167. 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));
  168. 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));
  169. 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));
  170. }
  171. }
  172. }
  173. //Manage output TProfiles
  174. TH2D *hPsi2EPRaw[NEtaGaps];
  175. TH2D *hPsi2EPRec[NEtaGaps];
  176. TH2D *hPsi2EPFlt[NEtaGaps];
  177. TH2D *hPsi3EPRaw[NEtaGaps];
  178. TH2D *hPsi3EPRec[NEtaGaps];
  179. TH2D *hPsi3EPFlt[NEtaGaps];
  180. for (int iGap=0;iGap<NEtaGaps;iGap++)
  181. {
  182. hPsi2EPRaw[iGap] = new TH2D(Form("hPsi2EPRaw%i",iGap),Form("hPsi2EPRaw%i",iGap),500,-5./2.,5./2.,9,-0.5,8.5);
  183. hPsi2EPRec[iGap] = new TH2D(Form("hPsi2EPRec%i",iGap),Form("hPsi2EPRec%i",iGap),500,-5./2.,5./2.,9,-0.5,8.5);
  184. hPsi2EPFlt[iGap] = new TH2D(Form("hPsi2EPFlt%i",iGap),Form("hPsi2EPFlt%i",iGap),500,-5./2.,5./2.,9,-0.5,8.5);
  185. hPsi3EPRaw[iGap] = new TH2D(Form("hPsi3EPRaw%i",iGap),Form("hPsi2EPRaw%i",iGap),500,-5./3.,5./3.,9,-0.5,8.5);
  186. hPsi3EPRec[iGap] = new TH2D(Form("hPsi3EPRec%i",iGap),Form("hPsi2EPRec%i",iGap),500,-5./3.,5./3.,9,-0.5,8.5);
  187. hPsi3EPFlt[iGap] = new TH2D(Form("hPsi3EPFlt%i",iGap),Form("hPsi2EPFlt%i",iGap),500,-5./3.,5./3.,9,-0.5,8.5);
  188. }
  189. //Counters
  190. Int_t fUsedTracks = 0;
  191. Int_t fQcountFull = 0;
  192. Int_t fQcountFullEast = 0;
  193. Int_t fQcountFullWest = 0;
  194. Int_t fQcountWest[NEtaGaps];
  195. Int_t fQcountEast[NEtaGaps];
  196. //Q-vectors
  197. TVector2 fQv2FullEP;
  198. TVector2 fQv2FullSP;
  199. TVector2 fQv3FullEP;
  200. TVector2 fQv3FullSP;
  201. TVector2 fQv2FullEPRaw;
  202. TVector2 fQv3FullEPRaw;
  203. TVector2 fQv2WestEP[NEtaGaps];
  204. TVector2 fQv2EastEP[NEtaGaps];
  205. TVector2 fQv2WestSP[NEtaGaps];
  206. TVector2 fQv2EastSP[NEtaGaps];
  207. TVector2 fQv3WestEP[NEtaGaps];
  208. TVector2 fQv3EastEP[NEtaGaps];
  209. TVector2 fQv3WestSP[NEtaGaps];
  210. TVector2 fQv3EastSP[NEtaGaps];
  211. TVector2 QvTMP;
  212. Double_t weight;
  213. Int_t VtxSign;
  214. Double_t psiShiftedTMP;
  215. Double_t PsiSin;
  216. Double_t PsiCos;
  217. Long64_t goodEventCounter = 0;
  218. Double_t Psi2EastEP[NEtaGaps], Psi3EastEP[NEtaGaps];
  219. Double_t Psi2WestEP[NEtaGaps], Psi3WestEP[NEtaGaps];
  220. Double_t Psi_ReCenter, dPsi;
  221. Int_t bin;
  222. Double_t res2EP[NEtaGaps], res3EP[NEtaGaps];
  223. // Loop over events
  224. for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
  225. {
  226. if ((iEvent+1) % 1000 == 0){
  227. std::cout << "Working on event #[" << (iEvent + 1)
  228. << "/" << events2read << "]" << std::endl;
  229. }
  230. Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
  231. if (!readEvent)
  232. {
  233. std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
  234. break;
  235. }
  236. // Retrieve femtoDst
  237. StFemtoDst *dst = femtoReader->femtoDst();
  238. // Retrieve event information
  239. StFemtoEvent *event = dst->event();
  240. if (!event)
  241. {
  242. std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
  243. break;
  244. }
  245. // Event selection
  246. if ( !isGoodEvent(event) ) continue;
  247. goodEventCounter++;
  248. TVector3 pVtx = event->primaryVertex();
  249. // Track analysis
  250. Int_t nTracks = dst->numberOfTracks();
  251. //Counters
  252. fUsedTracks = 0;
  253. fQcountFull = 0;
  254. fQcountFullEast = 0;
  255. fQcountFullWest = 0;
  256. psiShiftedTMP = 0.;
  257. PsiSin = 0.;
  258. PsiCos = 0.;
  259. fQv2FullEP.Set(0.,0.);
  260. fQv2FullEPRaw.Set(0.,0.);
  261. fQv2FullSP.Set(0.,0.);
  262. fQv3FullEP.Set(0.,0.);
  263. fQv3FullEPRaw.Set(0.,0.);
  264. fQv3FullSP.Set(0.,0.);
  265. QvTMP.Set(0.,0.);
  266. for (int i=0; i<NEtaGaps; i++)
  267. {
  268. fQcountWest[i] = 0;
  269. fQcountEast[i] = 0;
  270. fQv2WestEP[i].Set(0.,0.);
  271. fQv2EastEP[i].Set(0.,0.);
  272. fQv2WestSP[i].Set(0.,0.);
  273. fQv2EastSP[i].Set(0.,0.);
  274. fQv3WestEP[i].Set(0.,0.);
  275. fQv3EastEP[i].Set(0.,0.);
  276. fQv3WestSP[i].Set(0.,0.);
  277. fQv3EastSP[i].Set(0.,0.);
  278. }
  279. // Vertex sign
  280. VtxSign = (event->primaryVertex().Z() > 0) ? 0 : 1;
  281. // Track loop
  282. for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
  283. {
  284. // Retrieve i-th femto track
  285. StFemtoTrack *femtoTrack = dst->track(iTrk);
  286. if (!femtoTrack)
  287. continue;
  288. // Must be a primary track
  289. if (!femtoTrack->isPrimary())
  290. continue;
  291. //Track selection
  292. if (!isGoodTrack(femtoTrack, energy, pVtx)) continue;
  293. fUsedTracks++;
  294. //Fill recentered Q-vectors for all TPC (FULL)
  295. weight = GetWeight(femtoTrack);
  296. QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign]->GetBinContent(p_Qx_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  297. QvTMP.SetY(p_Qy_FULL_EP[0][VtxSign]->GetBinContent(p_Qy_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  298. fQv2FullEPRaw += weight*(CalcQvector(femtoTrack,2));
  299. fQv2FullEP += weight*(CalcQvector(femtoTrack,2) - QvTMP);
  300. QvTMP.Set(0.,0.);
  301. QvTMP.SetX(p_Qx_FULL_EP[1][VtxSign]->GetBinContent(p_Qx_FULL_EP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  302. QvTMP.SetY(p_Qy_FULL_EP[1][VtxSign]->GetBinContent(p_Qy_FULL_EP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  303. fQv3FullEPRaw += weight*(CalcQvector(femtoTrack,3));
  304. fQv3FullEP += weight*(CalcQvector(femtoTrack,3) - QvTMP);
  305. QvTMP.Set(0.,0.);
  306. QvTMP.SetX(p_Qx_FULL_SP[0][VtxSign]->GetBinContent(p_Qx_FULL_SP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  307. QvTMP.SetY(p_Qy_FULL_SP[0][VtxSign]->GetBinContent(p_Qy_FULL_SP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  308. fQv2FullSP += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
  309. QvTMP.Set(0.,0.);
  310. QvTMP.SetX(p_Qx_FULL_SP[1][VtxSign]->GetBinContent(p_Qx_FULL_SP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  311. QvTMP.SetY(p_Qy_FULL_SP[1][VtxSign]->GetBinContent(p_Qy_FULL_SP[1][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent9())));
  312. fQv3FullSP += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
  313. QvTMP.Set(0.,0.);
  314. fQcountFull++;
  315. if (femtoTrack->eta() >= 0.)
  316. {
  317. fQcountFullWest++;
  318. }
  319. else
  320. {
  321. fQcountFullEast++;
  322. }
  323. for (int iGap=0; iGap<NEtaGaps; iGap++)
  324. {
  325. //fill Qvectors from EAST arm
  326. if (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < cutEtaGap[iGap])
  327. {
  328. 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())));
  329. 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())));
  330. fQv2EastEP[iGap] += weight*(CalcQvector(femtoTrack,2) - QvTMP);
  331. QvTMP.Set(0.,0.);
  332. 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())));
  333. 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())));
  334. fQv3EastEP[iGap] += weight*(CalcQvector(femtoTrack,3) - QvTMP);
  335. QvTMP.Set(0.,0.);
  336. 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())));
  337. 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())));
  338. fQv2EastSP[iGap] += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
  339. QvTMP.Set(0.,0.);
  340. 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())));
  341. 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())));
  342. fQv2EastSP[iGap] += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
  343. QvTMP.Set(0.,0.);
  344. fQcountEast[iGap]++;
  345. }
  346. //fill Qvectors from WEST arm
  347. if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > cutEtaGap[iGap])
  348. {
  349. 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())));
  350. 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())));
  351. fQv2WestEP[iGap] += weight*(CalcQvector(femtoTrack,2) - QvTMP);
  352. QvTMP.Set(0.,0.);
  353. 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())));
  354. 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())));
  355. fQv3WestEP[iGap] += weight*(CalcQvector(femtoTrack,3) - QvTMP);
  356. QvTMP.Set(0.,0.);
  357. 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())));
  358. 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())));
  359. fQv2WestSP[iGap] += 1.*(CalcQvector(femtoTrack,2) - QvTMP);
  360. QvTMP.Set(0.,0.);
  361. 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())));
  362. 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())));
  363. fQv2WestSP[iGap] += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
  364. QvTMP.Set(0.,0.);
  365. fQcountWest[iGap]++;
  366. }
  367. }
  368. // Check if track has TOF signal
  369. if (femtoTrack->isTofTrack())
  370. {
  371. } //if( isTofTrack() )
  372. } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
  373. hPsi2EPRaw[0]->Fill((Double_t)TMath::ATan2(fQv2FullEPRaw.Y(),fQv2FullEPRaw.X())/2.,(Double_t)event->cent9());
  374. hPsi2EPRec[0]->Fill((Double_t)TMath::ATan2(fQv2FullEP.Y(),fQv2FullEP.X())/2.,(Double_t)event->cent9());
  375. hPsi3EPRaw[0]->Fill((Double_t)TMath::ATan2(fQv3FullEPRaw.Y(),fQv3FullEPRaw.X())/3.,(Double_t)event->cent9());
  376. hPsi3EPRec[0]->Fill((Double_t)TMath::ATan2(fQv3FullEP.Y(),fQv3FullEP.X())/3.,(Double_t)event->cent9());
  377. //Getting PsiEP with recentering + shift correction
  378. Psi_ReCenter=0.;
  379. Double_t mCos, mSin;
  380. Double_t Psi2EPFlat, Psi3EPFlat, Psi2EP, Psi3EP;
  381. Psi2EP = TMath::ATan2(fQv2FullEP.Y(),fQv2FullEP.X())/2.;
  382. Psi3EP = TMath::ATan2(fQv3FullEP.Y(),fQv3FullEP.X())/3.;
  383. dPsi = 0.;
  384. for (int iK=0; iK<NShiftOrderMax; iK++)
  385. {
  386. bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  387. mCos = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);
  388. bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  389. mSin = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);
  390. dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
  391. TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
  392. mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi2EP));
  393. }
  394. Psi2EPFlat = AngleShift(Psi2EP+dPsi,2.);
  395. dPsi = 0.;
  396. for (int iK=0; iK<NShiftOrderMax; iK++)
  397. {
  398. bin = p_Cos3_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  399. mCos = p_Cos3_FULL_EP[VtxSign][iK]->GetBinContent(bin);
  400. bin = p_Cos3_FULL_EP[VtxSign][iK]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  401. mSin = p_Cos3_FULL_EP[VtxSign][iK]->GetBinContent(bin);
  402. dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
  403. TMath::Cos(shiftOrderCoeff3.at(iK)*Psi3EP)+
  404. mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi3EP));
  405. }
  406. Psi3EPFlat = AngleShift(Psi3EP+dPsi,3.);
  407. hPsi2EPFlt[0]->Fill(Psi2EPFlat,(Double_t)event->cent9());
  408. hPsi3EPFlt[0]->Fill(Psi3EPFlat,(Double_t)event->cent9());
  409. for (int iGap=0; iGap<NEtaGaps; iGap++)
  410. {
  411. if (fQcountEast[iGap] > cutNcountQvGap &&
  412. fQcountWest[iGap] > cutNcountQvGap)
  413. {
  414. Psi2EastEP[iGap]=0.;
  415. Psi3EastEP[iGap]=0.;
  416. Psi2WestEP[iGap]=0.;
  417. Psi3WestEP[iGap]=0.;
  418. //EAST 2-nd order
  419. Psi_ReCenter = TMath::ATan2(fQv2EastEP[iGap].Y(),fQv2EastEP[iGap].X())/2.;
  420. dPsi = 0.;
  421. for (int iK=0; iK<NShiftOrderMax; iK++)
  422. {
  423. bin = p_Cos2_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  424. mCos = p_Cos2_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  425. bin = p_Cos2_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  426. mSin = p_Cos2_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  427. dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
  428. TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
  429. mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
  430. }
  431. Psi2EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
  432. //WEST 2-nd order
  433. Psi_ReCenter = TMath::ATan2(fQv2WestEP[iGap].Y(),fQv2WestEP[iGap].X())/2.;
  434. dPsi = 0.;
  435. for (int iK=0; iK<NShiftOrderMax; iK++)
  436. {
  437. bin = p_Cos2_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  438. mCos = p_Cos2_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  439. bin = p_Cos2_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  440. mSin = p_Cos2_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  441. dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
  442. TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
  443. mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
  444. }
  445. Psi2WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
  446. //EAST 3-nd order
  447. Psi_ReCenter = TMath::ATan2(fQv3EastEP[iGap].Y(),fQv3EastEP[iGap].X())/3.;
  448. dPsi = 0.;
  449. for (int iK=0; iK<NShiftOrderMax; iK++)
  450. {
  451. bin = p_Cos3_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  452. mCos = p_Cos3_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  453. bin = p_Cos3_EAST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  454. mSin = p_Cos3_EAST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  455. dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
  456. TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
  457. mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
  458. }
  459. Psi3EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
  460. //WEST 3-nd order
  461. Psi_ReCenter = TMath::ATan2(fQv3WestEP[iGap].Y(),fQv3WestEP[iGap].X())/3.;
  462. dPsi = 0.;
  463. for (int iK=0; iK<NShiftOrderMax; iK++)
  464. {
  465. bin = p_Cos3_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  466. mCos = p_Cos3_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  467. bin = p_Cos3_WEST_EP[VtxSign][iK][iGap]->FindBin((Double_t) event->runId(),(Double_t) event->cent9());
  468. mSin = p_Cos3_WEST_EP[VtxSign][iK][iGap]->GetBinContent(bin);
  469. dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
  470. TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
  471. mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
  472. }
  473. Psi3WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
  474. }
  475. }
  476. } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
  477. femtoReader->Finish();
  478. TFile *output = new TFile(outFile,"recreate");
  479. for (int iGap=0;iGap<NEtaGaps;iGap++)
  480. {
  481. hPsi2EPRaw[iGap]->Write();
  482. hPsi2EPRec[iGap]->Write();
  483. hPsi2EPFlt[iGap]->Write();
  484. hPsi3EPRaw[iGap]->Write();
  485. hPsi3EPRec[iGap]->Write();
  486. hPsi3EPFlt[iGap]->Write();
  487. }
  488. output->Close();
  489. std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
  490. << std::endl;
  491. std::cout << "Good events: " << goodEventCounter << std::endl;
  492. timer.Stop();
  493. timer.Print();
  494. }
  495. Bool_t isGoodEvent(StFemtoEvent *const &event)
  496. {
  497. if (!event) return false;
  498. if (event == nullptr) return false;
  499. if (event->primaryVertex().Perp() > cutVtxR) return false;
  500. if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy)) return false;
  501. if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz) return false;
  502. return true;
  503. }
  504. Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
  505. {
  506. if (!track) return false;
  507. // if (!track->isPrimary()) return false;
  508. if (track->nHitsFit() < cutNhits) return false;
  509. if (track->nHitsPoss() <= cutNhitsPoss) return false;
  510. if ((Double_t) track->nHitsFit()/track->nHitsPoss() < cutNhitsRatio) return false;
  511. if (TMath::Abs(track->eta()) >= cutEta) return false;
  512. if (track->pt() <= cutPtMin.at(_energy)) return false;
  513. if (track->pt() > cutPtMax) return false;
  514. if (track->ptot() > cutPMax) return false;
  515. if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
  516. return true;
  517. }
  518. Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
  519. {
  520. if (!track) return false;
  521. // if (!track->isPrimary()) return false;
  522. if (track->nHitsFit() < cutNhits) return false;
  523. if (track->nHitsPoss() <= cutNhitsPoss) return false;
  524. if ((Double_t) track->nHitsFit()/track->nHitsPoss() < cutNhitsRatio) return false;
  525. if (TMath::Abs(track->eta()) >= cutEta) return false;
  526. if (track->pt() <= cutPtMin.at(_energy)) return false;
  527. //if (track->pt() > cutPtMax) return false;
  528. if (track->ptot() > cutPMax) return false;
  529. if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
  530. return true;
  531. }
  532. Double_t GetWeight(StFemtoTrack *const &track)
  533. {
  534. Double_t weight;
  535. if (track->pt() <= cutPtWeightEP)
  536. {
  537. weight = track->pt();
  538. }
  539. else
  540. {
  541. weight = cutPtWeightEP;
  542. }
  543. return weight;
  544. }
  545. TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
  546. {
  547. TVector2 qv(0.,0.);
  548. qv.Set(TMath::Cos(_harm*track->phi()),TMath::Sin(_harm*track->phi()));
  549. return qv;
  550. }
  551. Double_t AngleShift(Double_t Psi, Double_t order)
  552. {
  553. Double_t PsiCorr = Psi;
  554. if (PsiCorr > TMath::Pi()/order)
  555. {
  556. PsiCorr = Psi - 2.*TMath::Pi()/order;
  557. }
  558. if (PsiCorr < -1.*TMath::Pi()/order)
  559. {
  560. PsiCorr = Psi + 2.*TMath::Pi()/order;
  561. }
  562. return PsiCorr;
  563. }