FemtoDstAnalyzer_ZDC_ShiftCorr.C 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375
  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_ShiftCorr(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
  50. const Char_t *outFile = "./oReCenteringTest.root",
  51. const Char_t *recZDC = "")
  52. {
  53. std::cout << "Hi! Lets do some physics, Master!" << std::endl;
  54. TStopwatch timer;
  55. timer.Start();
  56. #if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
  57. gSystem->Load("../libStFemtoDst.so");
  58. #endif
  59. //Manage recentering
  60. TProfile2D *p_ZDC_Qx_Full_EP[fNharmonics+1][fArraySize];
  61. TProfile2D *p_ZDC_Qy_Full_EP[fNharmonics+1][fArraySize];
  62. TProfile2D *p_ZDC_Qx_East_EP[fNharmonics+1][fArraySize];
  63. TProfile2D *p_ZDC_Qy_East_EP[fNharmonics+1][fArraySize];
  64. TProfile2D *p_ZDC_Qx_West_EP[fNharmonics+1][fArraySize];
  65. TProfile2D *p_ZDC_Qy_West_EP[fNharmonics+1][fArraySize];
  66. TFile *fiReCentZDC = new TFile(recZDC, "read");
  67. fiReCentZDC->cd();
  68. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  69. {
  70. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  71. {
  72. p_ZDC_Qx_Full_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
  73. p_ZDC_Qy_Full_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
  74. p_ZDC_Qx_East_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
  75. p_ZDC_Qy_East_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
  76. p_ZDC_Qx_West_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
  77. p_ZDC_Qy_West_EP[iHarm][iVtx] = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
  78. }
  79. }
  80. //Manage output
  81. TProfile2D *p_ZDC_Cos_Full_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  82. TProfile2D *p_ZDC_Sin_Full_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  83. TProfile2D *p_ZDC_Cos_East_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  84. TProfile2D *p_ZDC_Sin_East_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  85. TProfile2D *p_ZDC_Cos_West_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  86. TProfile2D *p_ZDC_Sin_West_EP[fNharmonics+1][fArraySize][NShiftOrderMax*3];
  87. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  88. {
  89. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  90. {
  91. for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
  92. {
  93. p_ZDC_Cos_Full_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  94. p_ZDC_Sin_Full_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  95. p_ZDC_Cos_East_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  96. p_ZDC_Sin_East_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  97. p_ZDC_Cos_West_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  98. p_ZDC_Sin_West_EP[iHarm][iVtx][iShift] = new TProfile2D(Form("p_ZDC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), Form("p_ZDC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx), fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
  99. }
  100. }
  101. }
  102. StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
  103. femtoReader->Init();
  104. // This is a way if you want to spead up IO
  105. std::cout << "Explicit read status for some branches" << std::endl;
  106. femtoReader->SetStatus("*", 0);
  107. femtoReader->SetStatus("Event", 1);
  108. femtoReader->SetStatus("Track", 1);
  109. std::cout << "Status has been set" << std::endl;
  110. std::cout << "Now I know what to read, Master!" << std::endl;
  111. if (!femtoReader->chain())
  112. {
  113. std::cout << "No chain has been found." << std::endl;
  114. }
  115. Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  116. std::cout << "eventsInTree: " << eventsInTree << std::endl;
  117. Long64_t events2read = femtoReader->chain()->GetEntries();
  118. std::cout << "Number of events to read: " << events2read << std::endl;
  119. Int_t VtxSign;
  120. // Loop over events
  121. for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
  122. {
  123. if ((iEvent + 1) % 1000 == 0)
  124. {
  125. std::cout << "Working on event #[" << (iEvent + 1)
  126. << "/" << events2read << "]" << std::endl;
  127. }
  128. Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
  129. if (!readEvent)
  130. {
  131. std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
  132. break;
  133. }
  134. // Retrieve femtoDst
  135. StFemtoDst *dst = femtoReader->femtoDst();
  136. // Retrieve event information
  137. StFemtoEvent *event = dst->event();
  138. if (!event)
  139. {
  140. std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
  141. break;
  142. }
  143. // Event selection
  144. if (!isGoodEvent(event))
  145. continue;
  146. VtxSign = GetVzBin(event->primaryVertex().Z());
  147. if (VtxSign == -1)
  148. continue;
  149. if (event->zdcSumAdcEast() < 1.) continue;
  150. if (event->zdcSumAdcWest() < 1.) continue;
  151. Double_t qxEast, qyEast, adcxEast, adcyEast;
  152. Double_t qxWest, qyWest, adcxWest, adcyWest;
  153. Double_t qxFull, qyFull, adcxFull, adcyFull;
  154. qxEast = 0.;
  155. qyEast = 0.;
  156. adcxEast = 0.;
  157. adcyEast = 0.;
  158. qxWest = 0.;
  159. qyWest = 0.;
  160. adcxWest = 0.;
  161. adcyWest = 0.;
  162. qxFull = 0.;
  163. qyFull = 0.;
  164. adcxFull = 0.;
  165. adcyFull = 0.;
  166. Double_t qxRecEast, qyRecEast;
  167. Double_t qxRecWest, qyRecWest;
  168. Double_t qxRecFull, qyRecFull;
  169. // Vertical ZDC-SMD
  170. for (int iStrip=0; iStrip<fNZdcSmdStripsVertical; iStrip++)
  171. {
  172. qxEast += event->zdcSmdEastVertical(iStrip) * GetZDCPosition(0,0,iStrip);
  173. adcxEast += event->zdcSmdEastVertical(iStrip);
  174. qxWest += event->zdcSmdWestVertical(iStrip) * GetZDCPosition(1,0,iStrip);
  175. adcxWest += event->zdcSmdWestVertical(iStrip);
  176. }
  177. // Horizontal ZDC-SMD
  178. for (int iStrip=0; iStrip<fNZdcSmdStripsHorizontal; iStrip++)
  179. {
  180. qyEast += event->zdcSmdEastHorizontal(iStrip) * GetZDCPosition(0,1,iStrip);
  181. adcyEast += event->zdcSmdEastHorizontal(iStrip);
  182. qyWest += event->zdcSmdWestHorizontal(iStrip) * GetZDCPosition(1,1,iStrip);
  183. adcyWest += event->zdcSmdWestHorizontal(iStrip);
  184. }
  185. if (adcxEast > 0) qxEast /= adcxEast;
  186. else qxEast = -999.;
  187. if (adcyEast > 0) qyEast /= adcyEast;
  188. else qyEast = -999.;
  189. if (adcxWest > 0) qxWest /= adcxWest;
  190. else qxWest = -999.;
  191. if (adcyWest > 0) qyWest /= adcyWest;
  192. else qyWest = -999.;
  193. // Full ZDC-SMD Q-vector
  194. if (qxEast == -999. || qxWest == -999.) qxFull = -999.;
  195. else qxFull = qxEast - qxWest;
  196. if (qyEast == -999. || qyWest == -999.) qyFull = -999.;
  197. else qyFull = qyEast - qyWest;
  198. if (qxEast != -999.) qxRecEast = qxEast - (Double_t)p_ZDC_Qx_East_EP[0][VtxSign]->GetBinContent(p_ZDC_Qx_East_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
  199. if (qyEast != -999.) qyRecEast = qyEast - (Double_t)p_ZDC_Qy_East_EP[0][VtxSign]->GetBinContent(p_ZDC_Qy_East_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
  200. if (qxWest != -999.) qxRecWest = qxWest - (Double_t)p_ZDC_Qx_West_EP[0][VtxSign]->GetBinContent(p_ZDC_Qx_West_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
  201. if (qyWest != -999.) qyRecWest = qyWest - (Double_t)p_ZDC_Qy_West_EP[0][VtxSign]->GetBinContent(p_ZDC_Qy_West_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
  202. if (qxFull != -999.) qxRecFull = qxFull - (Double_t)p_ZDC_Qx_Full_EP[0][VtxSign]->GetBinContent(p_ZDC_Qx_Full_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
  203. if (qyFull != -999.) qyRecFull = qyFull - (Double_t)p_ZDC_Qy_Full_EP[0][VtxSign]->GetBinContent(p_ZDC_Qy_Full_EP[0][VtxSign]->FindBin(event->runId(), event->cent16()));
  204. if (qxEast == -999.) qxRecEast = qxEast;
  205. if (qyEast == -999.) qyRecEast = qyEast;
  206. if (qxWest == -999.) qxRecWest = qxWest;
  207. if (qyWest == -999.) qyRecWest = qyWest;
  208. if (qxFull == -999.) qxRecFull = qxFull;
  209. if (qyFull == -999.) qyRecFull = qyFull;
  210. Double_t PsiSin = 0., PsiCos = 0., psiShiftedTMP = 0.;
  211. TVector2 qRecEast[fNharmonics+1], qRecWest[fNharmonics+1], qRecFull[fNharmonics+1];
  212. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  213. {
  214. qRecEast[iHarm].Set(qxRecEast,qyRecEast);
  215. qRecWest[iHarm].Set(qxRecWest,qyRecWest);
  216. qRecFull[iHarm].Set(qxRecFull,qyRecFull);
  217. for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
  218. {
  219. //Full ZDC
  220. if (qRecFull[iHarm].Mod())
  221. {
  222. psiShiftedTMP = qRecFull[iHarm].Phi() / (iHarm + 1.);
  223. if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
  224. }
  225. PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
  226. PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
  227. p_ZDC_Cos_Full_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
  228. p_ZDC_Sin_Full_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);
  229. //East ZDC
  230. if (qRecEast[iHarm].Mod())
  231. {
  232. psiShiftedTMP = qRecEast[iHarm].Phi() / (iHarm + 1.);
  233. if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
  234. }
  235. PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
  236. PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
  237. p_ZDC_Cos_East_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
  238. p_ZDC_Sin_East_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);
  239. //West ZDC
  240. if (qRecWest[iHarm].Mod())
  241. {
  242. psiShiftedTMP = qRecWest[iHarm].Phi() / (iHarm + 1.);
  243. if (psiShiftedTMP < 0.0) {psiShiftedTMP += 2.*TMath::Pi()/(iHarm+1.);}
  244. }
  245. PsiSin = TMath::Sin((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
  246. PsiCos = TMath::Cos((iHarm+1.)*(iShift+1.) * psiShiftedTMP);
  247. p_ZDC_Cos_West_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiCos);
  248. p_ZDC_Sin_West_EP[iHarm][VtxSign][iShift]->Fill(event->runId(), event->cent16(), PsiSin);
  249. }
  250. }
  251. } // end of the event loop
  252. femtoReader->Finish();
  253. TFile *output = new TFile(outFile, "recreate");
  254. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  255. {
  256. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  257. {
  258. for (int iShift = 0; iShift < NShiftOrderMax*3; iShift++)
  259. {
  260. p_ZDC_Cos_Full_EP[iHarm][iVtx][iShift]->Write();
  261. p_ZDC_Sin_Full_EP[iHarm][iVtx][iShift]->Write();
  262. p_ZDC_Cos_East_EP[iHarm][iVtx][iShift]->Write();
  263. p_ZDC_Sin_East_EP[iHarm][iVtx][iShift]->Write();
  264. p_ZDC_Cos_West_EP[iHarm][iVtx][iShift]->Write();
  265. p_ZDC_Sin_West_EP[iHarm][iVtx][iShift]->Write();
  266. }
  267. }
  268. }
  269. output->Close();
  270. std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
  271. << std::endl;
  272. timer.Stop();
  273. timer.Print();
  274. }
  275. Bool_t isGoodEvent(StFemtoEvent *const &event)
  276. {
  277. if (!event)
  278. return false;
  279. if (event == nullptr)
  280. return false;
  281. if (event->primaryVertex().Perp() > cutVtxR)
  282. return false;
  283. if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy))
  284. return false;
  285. if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz)
  286. return false;
  287. return true;
  288. }
  289. Int_t GetVzBin(Double_t vtxZ)
  290. {
  291. for (Int_t i = 0; i < fArraySize; i++)
  292. {
  293. if (vtxZ >= VtxBins[i] && vtxZ < VtxBins[i + 1])
  294. return i;
  295. }
  296. return -1;
  297. }
  298. Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip)
  299. // Get position of each slat;strip starts from 0
  300. {
  301. Double_t zdcsmd_x[7] = {0.5,2,3.5,5,6.5,8,9.5};
  302. Double_t zdcsmd_y[8] = {1.25,3.25,5.25,7.25,9.25,11.25,13.25,15.25};
  303. Double_t mZDCSMDCenterex = 4.72466;
  304. Double_t mZDCSMDCenterey = 5.53629;
  305. Double_t mZDCSMDCenterwx = 4.39604;
  306. Double_t mZDCSMDCenterwy = 5.19968;
  307. if(eastwest==0 && verthori==0) return zdcsmd_x[strip]-mZDCSMDCenterex;
  308. if(eastwest==1 && verthori==0) return mZDCSMDCenterwx-zdcsmd_x[strip];
  309. if(eastwest==0 && verthori==1) return zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterey;
  310. if(eastwest==1 && verthori==1) return zdcsmd_y[strip]/sqrt(2.)-mZDCSMDCenterwy;
  311. return -999.;
  312. }