FemtoDstAnalyzer_ZDC_ReCentering.C 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288
  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. #include "TRandom3.h"
  31. // FemtoDst headers
  32. #include "StFemtoDstReader.h"
  33. #include "StFemtoDst.h"
  34. #include "StFemtoEvent.h"
  35. #include "StFemtoTrack.h"
  36. // Load libraries (for ROOT_VERSTION_CODE >= 393215)
  37. #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
  38. R__LOAD_LIBRARY(/ mnt / pool / rhic / 4 / parfenovpeter / STAR / build / libStFemtoDst.so)
  39. #endif
  40. #include "Constants.h"
  41. // inFile - is a name of name.FemtoDst.root file or a name
  42. // of a name.lis(t) files, that contains a list of
  43. // name1.FemtoDst.root, name2.FemtoDst.root, ... files
  44. //Used functions (see them below)
  45. Bool_t isGoodEvent(StFemtoEvent *const &event);
  46. Int_t GetVzBin(Double_t vtxZ);
  47. Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip);
  48. //_________________
  49. void FemtoDstAnalyzer_ZDC_ReCentering(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
  50. const Char_t *outFile = "./oReCenteringTest.root")
  51. {
  52. std::cout << "Hi! Lets do some physics, Master!" << std::endl;
  53. TStopwatch timer;
  54. timer.Start();
  55. #if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
  56. gSystem->Load("../libStFemtoDst.so");
  57. #endif
  58. //Manage output
  59. TProfile2D *p_ZDC_Qx_Full_EP[fNharmonics+1][fArraySize];
  60. TProfile2D *p_ZDC_Qy_Full_EP[fNharmonics+1][fArraySize];
  61. TProfile2D *p_ZDC_Qx_East_EP[fNharmonics+1][fArraySize];
  62. TProfile2D *p_ZDC_Qy_East_EP[fNharmonics+1][fArraySize];
  63. TProfile2D *p_ZDC_Qx_West_EP[fNharmonics+1][fArraySize];
  64. TProfile2D *p_ZDC_Qy_West_EP[fNharmonics+1][fArraySize];
  65. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  66. {
  67. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  68. {
  69. p_ZDC_Qx_Full_EP[iHarm][iVtx] = new TProfile2D(Form("p_ZDC_Qx_Full_EP%i_vtx%i", iHarm, iVtx), Form("p_ZDC_Qx_Full_EP%i_vtx%i", iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  70. p_ZDC_Qy_Full_EP[iHarm][iVtx] = new TProfile2D(Form("p_ZDC_Qy_Full_EP%i_vtx%i", iHarm, iVtx), Form("p_ZDC_Qy_Full_EP%i_vtx%i", iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  71. p_ZDC_Qx_East_EP[iHarm][iVtx] = new TProfile2D(Form("p_ZDC_Qx_East_EP%i_vtx%i", iHarm, iVtx), Form("p_ZDC_Qx_East_EP%i_vtx%i", iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  72. p_ZDC_Qy_East_EP[iHarm][iVtx] = new TProfile2D(Form("p_ZDC_Qy_East_EP%i_vtx%i", iHarm, iVtx), Form("p_ZDC_Qy_East_EP%i_vtx%i", iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  73. p_ZDC_Qx_West_EP[iHarm][iVtx] = new TProfile2D(Form("p_ZDC_Qx_West_EP%i_vtx%i", iHarm, iVtx), Form("p_ZDC_Qx_West_EP%i_vtx%i", iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  74. p_ZDC_Qy_West_EP[iHarm][iVtx] = new TProfile2D(Form("p_ZDC_Qy_West_EP%i_vtx%i", iHarm, iVtx), Form("p_ZDC_Qy_West_EP%i_vtx%i", iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  75. }
  76. }
  77. StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
  78. femtoReader->Init();
  79. // This is a way if you want to spead up IO
  80. std::cout << "Explicit read status for some branches" << std::endl;
  81. femtoReader->SetStatus("*", 0);
  82. femtoReader->SetStatus("Event", 1);
  83. femtoReader->SetStatus("Track", 1);
  84. std::cout << "Status has been set" << std::endl;
  85. std::cout << "Now I know what to read, Master!" << std::endl;
  86. if (!femtoReader->chain())
  87. {
  88. std::cout << "No chain has been found." << std::endl;
  89. }
  90. Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  91. std::cout << "eventsInTree: " << eventsInTree << std::endl;
  92. Long64_t events2read = femtoReader->chain()->GetEntries();
  93. std::cout << "Number of events to read: " << events2read << std::endl;
  94. Int_t VtxSign;
  95. // Loop over events
  96. for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
  97. {
  98. if ((iEvent + 1) % 1000 == 0)
  99. {
  100. std::cout << "Working on event #[" << (iEvent + 1)
  101. << "/" << events2read << "]" << std::endl;
  102. }
  103. Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
  104. if (!readEvent)
  105. {
  106. std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
  107. break;
  108. }
  109. // Retrieve femtoDst
  110. StFemtoDst *dst = femtoReader->femtoDst();
  111. // Retrieve event information
  112. StFemtoEvent *event = dst->event();
  113. if (!event)
  114. {
  115. std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
  116. break;
  117. }
  118. // Event selection
  119. if (!isGoodEvent(event))
  120. continue;
  121. VtxSign = GetVzBin(event->primaryVertex().Z());
  122. if (VtxSign == -1)
  123. continue;
  124. if (event->zdcSumAdcEast() < 1.) continue;
  125. if (event->zdcSumAdcWest() < 1.) continue;
  126. Double_t qxEast, qyEast, adcxEast, adcyEast;
  127. Double_t qxWest, qyWest, adcxWest, adcyWest;
  128. Double_t qxFull, qyFull, adcxFull, adcyFull;
  129. qxEast = 0.;
  130. qyEast = 0.;
  131. adcxEast = 0.;
  132. adcyEast = 0.;
  133. qxWest = 0.;
  134. qyWest = 0.;
  135. adcxWest = 0.;
  136. adcyWest = 0.;
  137. qxFull = 0.;
  138. qyFull = 0.;
  139. adcxFull = 0.;
  140. adcyFull = 0.;
  141. // Vertical ZDC-SMD
  142. for (int iStrip=0; iStrip<fNZdcSmdStripsVertical; iStrip++)
  143. {
  144. qxEast += event->zdcSmdEastVertical(iStrip) * GetZDCPosition(0,0,iStrip);
  145. adcxEast += event->zdcSmdEastVertical(iStrip);
  146. qxWest += event->zdcSmdWestVertical(iStrip) * GetZDCPosition(1,0,iStrip);
  147. adcxWest += event->zdcSmdWestVertical(iStrip);
  148. }
  149. // Horizontal ZDC-SMD
  150. for (int iStrip=0; iStrip<fNZdcSmdStripsHorizontal; iStrip++)
  151. {
  152. qyEast += event->zdcSmdEastHorizontal(iStrip) * GetZDCPosition(0,1,iStrip);
  153. adcyEast += event->zdcSmdEastHorizontal(iStrip);
  154. qyWest += event->zdcSmdWestHorizontal(iStrip) * GetZDCPosition(1,1,iStrip);
  155. adcyWest += event->zdcSmdWestHorizontal(iStrip);
  156. }
  157. if (adcxEast > 0) qxEast /= adcxEast;
  158. else qxEast = -999.;
  159. if (adcyEast > 0) qyEast /= adcyEast;
  160. else qyEast = -999.;
  161. if (adcxWest > 0) qxWest /= adcxWest;
  162. else qxWest = -999.;
  163. if (adcyWest > 0) qyWest /= adcyWest;
  164. else qyWest = -999.;
  165. // Full ZDC-SMD Q-vector
  166. if (qxEast == -999. || qxWest == -999.) qxFull = -999.;
  167. else qxFull = qxEast - qxWest;
  168. if (qyEast == -999. || qyWest == -999.) qyFull = -999.;
  169. else qyFull = qyEast - qyWest;
  170. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  171. {
  172. p_ZDC_Qx_Full_EP[iHarm][VtxSign]->Fill(event->runId(), event->cent16(), qxFull);
  173. p_ZDC_Qy_Full_EP[iHarm][VtxSign]->Fill(event->runId(), event->cent16(), qyFull);
  174. p_ZDC_Qx_East_EP[iHarm][VtxSign]->Fill(event->runId(), event->cent16(), qxEast);
  175. p_ZDC_Qy_East_EP[iHarm][VtxSign]->Fill(event->runId(), event->cent16(), qyEast);
  176. p_ZDC_Qx_West_EP[iHarm][VtxSign]->Fill(event->runId(), event->cent16(), qxWest);
  177. p_ZDC_Qy_West_EP[iHarm][VtxSign]->Fill(event->runId(), event->cent16(), qyWest);
  178. }
  179. } // end of the event loop
  180. femtoReader->Finish();
  181. TFile *output = new TFile(outFile, "recreate");
  182. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  183. {
  184. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  185. {
  186. p_ZDC_Qx_Full_EP[iHarm][iVtx]->Write();
  187. p_ZDC_Qy_Full_EP[iHarm][iVtx]->Write();
  188. p_ZDC_Qx_East_EP[iHarm][iVtx]->Write();
  189. p_ZDC_Qy_East_EP[iHarm][iVtx]->Write();
  190. p_ZDC_Qx_West_EP[iHarm][iVtx]->Write();
  191. p_ZDC_Qy_West_EP[iHarm][iVtx]->Write();
  192. }
  193. }
  194. output->Close();
  195. std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
  196. << std::endl;
  197. timer.Stop();
  198. timer.Print();
  199. }
  200. Bool_t isGoodEvent(StFemtoEvent *const &event)
  201. {
  202. if (!event)
  203. return false;
  204. if (event == nullptr)
  205. return false;
  206. if (event->primaryVertex().Perp() > cutVtxR)
  207. return false;
  208. if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy))
  209. return false;
  210. if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz)
  211. return false;
  212. return true;
  213. }
  214. Int_t GetVzBin(Double_t vtxZ)
  215. {
  216. for (Int_t i = 0; i < fArraySize; i++)
  217. {
  218. if (vtxZ >= VtxBins[i] && vtxZ < VtxBins[i + 1])
  219. return i;
  220. }
  221. return -1;
  222. }
  223. Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip)
  224. // Get position of each slat;strip starts from 0
  225. {
  226. Double_t zdcsmd_x[7] = {0.5,2,3.5,5,6.5,8,9.5};
  227. Double_t zdcsmd_y[8] = {1.25,3.25,5.25,7.25,9.25,11.25,13.25,15.25};
  228. Double_t mZDCSMDCenterex = 4.72466;
  229. Double_t mZDCSMDCenterey = 5.53629;
  230. Double_t mZDCSMDCenterwx = 4.39604;
  231. Double_t mZDCSMDCenterwy = 5.19968;
  232. if(eastwest==0 && verthori==0) return zdcsmd_x[strip]-mZDCSMDCenterex;
  233. if(eastwest==1 && verthori==0) return mZDCSMDCenterwx-zdcsmd_x[strip];
  234. if(eastwest==0 && verthori==1) return zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterey;
  235. if(eastwest==1 && verthori==1) return zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterwy;
  236. return -999.;
  237. }