FemtoDstAnalyzer_FlowPIDHadrons.C 101 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050105110521053105410551056105710581059106010611062106310641065106610671068106910701071107210731074107510761077107810791080108110821083108410851086108710881089109010911092109310941095109610971098109911001101110211031104110511061107110811091110111111121113111411151116111711181119112011211122112311241125112611271128112911301131113211331134113511361137113811391140114111421143114411451146114711481149115011511152115311541155115611571158115911601161116211631164116511661167116811691170117111721173117411751176117711781179118011811182118311841185118611871188118911901191119211931194119511961197119811991200120112021203120412051206120712081209121012111212121312141215121612171218121912201221122212231224122512261227122812291230123112321233123412351236123712381239124012411242124312441245124612471248124912501251125212531254125512561257125812591260126112621263126412651266126712681269127012711272127312741275127612771278127912801281128212831284128512861287128812891290129112921293129412951296129712981299130013011302130313041305130613071308130913101311131213131314131513161317131813191320132113221323132413251326132713281329133013311332133313341335133613371338133913401341134213431344134513461347134813491350135113521353135413551356135713581359136013611362136313641365136613671368136913701371137213731374137513761377137813791380138113821383138413851386138713881389139013911392139313941395139613971398139914001401140214031404140514061407140814091410141114121413141414151416141714181419142014211422142314241425142614271428142914301431143214331434143514361437143814391440144114421443144414451446144714481449145014511452145314541455145614571458145914601461146214631464146514661467146814691470147114721473147414751476147714781479148014811482148314841485148614871488148914901491149214931494149514961497149814991500150115021503150415051506150715081509151015111512151315141515151615171518151915201521152215231524152515261527152815291530153115321533153415351536153715381539154015411542154315441545154615471548154915501551155215531554155515561557155815591560156115621563156415651566156715681569157015711572157315741575157615771578157915801581158215831584158515861587158815891590159115921593159415951596159715981599160016011602160316041605160616071608160916101611161216131614161516161617161816191620162116221623162416251626162716281629163016311632163316341635163616371638163916401641164216431644164516461647164816491650165116521653165416551656165716581659166016611662166316641665166616671668166916701671167216731674167516761677167816791680168116821683168416851686168716881689169016911692169316941695169616971698169917001701170217031704170517061707170817091710171117121713171417151716171717181719172017211722172317241725172617271728172917301731173217331734173517361737173817391740174117421743174417451746174717481749175017511752175317541755175617571758175917601761176217631764176517661767176817691770177117721773177417751776177717781779178017811782178317841785178617871788178917901791179217931794179517961797179817991800180118021803180418051806180718081809181018111812181318141815181618171818181918201821182218231824182518261827182818291830183118321833183418351836183718381839184018411842184318441845184618471848184918501851185218531854185518561857185818591860186118621863186418651866186718681869187018711872187318741875187618771878187918801881188218831884188518861887188818891890189118921893189418951896189718981899190019011902190319041905190619071908190919101911191219131914191519161917191819191920192119221923192419251926192719281929193019311932193319341935193619371938193919401941194219431944194519461947194819491950195119521953195419551956195719581959196019611962196319641965196619671968196919701971197219731974197519761977197819791980198119821983198419851986198719881989199019911992199319941995199619971998199920002001200220032004200520062007200820092010201120122013201420152016201720182019202020212022202320242025202620272028202920302031203220332034203520362037203820392040204120422043204420452046204720482049205020512052205320542055205620572058205920602061206220632064206520662067206820692070207120722073207420752076207720782079208020812082208320842085208620872088208920902091209220932094209520962097209820992100210121022103210421052106210721082109211021112112211321142115211621172118211921202121212221232124212521262127212821292130213121322133213421352136213721382139214021412142214321442145214621472148214921502151215221532154215521562157215821592160216121622163216421652166216721682169217021712172217321742175217621772178217921802181218221832184218521862187218821892190219121922193219421952196219721982199220022012202220322042205220622072208220922102211221222132214221522162217221822192220222122222223222422252226222722282229223022312232223322342235223622372238223922402241224222432244224522462247224822492250225122522253225422552256225722582259226022612262226322642265226622672268226922702271227222732274227522762277227822792280228122822283228422852286228722882289229022912292229322942295229622972298229923002301230223032304230523062307230823092310231123122313231423152316231723182319232023212322232323242325232623272328232923302331233223332334233523362337233823392340234123422343234423452346234723482349235023512352235323542355235623572358235923602361236223632364236523662367236823692370237123722373237423752376237723782379238023812382238323842385238623872388238923902391239223932394239523962397239823992400240124022403240424052406240724082409241024112412241324142415241624172418241924202421242224232424242524262427242824292430243124322433243424352436243724382439244024412442244324442445244624472448244924502451245224532454245524562457245824592460246124622463246424652466246724682469247024712472247324742475247624772478247924802481248224832484248524862487248824892490249124922493249424952496249724982499250025012502250325042505250625072508250925102511251225132514251525162517251825192520252125222523252425252526252725282529253025312532253325342535253625372538253925402541254225432544254525462547254825492550255125522553255425552556255725582559256025612562256325642565256625672568256925702571257225732574257525762577257825792580258125822583258425852586258725882589259025912592259325942595259625972598259926002601260226032604260526062607260826092610261126122613261426152616261726182619262026212622262326242625262626272628262926302631
  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 "TF1.h"
  28. #include "TMath.h"
  29. #include "TProfile2D.h"
  30. #include "TStopwatch.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. Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
  47. Bool_t isGoodTrackFlowPID(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
  48. Bool_t isGoodTrackFlowPIDsys_Nhits(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount);
  49. Bool_t isGoodTrackFlowPIDsys_DCA(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount);
  50. Bool_t isGoodTrackFlowPIDsys_Eta(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount);
  51. Bool_t isGoodPID(StFemtoTrack *const &track);
  52. Double_t GetWeight(StFemtoTrack *const &track);
  53. TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
  54. Double_t AngleShift(Double_t Psi, Double_t order);
  55. Int_t GetVzBin(Double_t vtxZ);
  56. //Int_t DoPID(StFemtoTrack *const &track);
  57. Double_t GetSigmWidth(StFemtoTrack *const &track, Int_t PID);
  58. Double_t GetMeanNsig(StFemtoTrack *const &track, const Int_t PID);
  59. Double_t GetWidthM2(StFemtoTrack *const &track, const Int_t PID);
  60. Double_t GetRotAngle(StFemtoTrack *const &track);
  61. std::pair<Double_t, Double_t> GetNewXY(StFemtoTrack *const &track);
  62. Bool_t isMesonXY(Double_t x, Double_t y);
  63. Double_t GetSigmNewX(StFemtoTrack *const &track, Int_t PID);
  64. Double_t GetSigmNewY(StFemtoTrack *const &track, Int_t PID);
  65. Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);
  66. Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip);
  67. //_________________
  68. void FemtoDstAnalyzer_FlowPIDHadrons(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
  69. const Char_t *outFile = "./oFlowPIDTest.root",
  70. const Char_t *reCentFile = "",
  71. const Char_t *shiftFile = "",
  72. const Char_t *resFile = "",
  73. const Char_t *gainBBC = "",
  74. const Char_t *reCentFileBBC = "",
  75. const Char_t *shiftFileBBC = "",
  76. const Char_t *reCentFileZDC = "",
  77. const Char_t *shiftFileZDC = "",
  78. const Char_t *resFitFile = "",
  79. const Char_t *pidFile = "")
  80. {
  81. std::cout << "Hi! Lets do some physics, Master!" << std::endl;
  82. TStopwatch timer;
  83. timer.Start();
  84. #if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
  85. gSystem->Load("../libStFemtoDst.so");
  86. #endif
  87. StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
  88. femtoReader->Init();
  89. // This is a way if you want to spead up IO
  90. std::cout << "Explicit read status for some branches" << std::endl;
  91. femtoReader->SetStatus("*", 0);
  92. femtoReader->SetStatus("Event", 1);
  93. femtoReader->SetStatus("Track", 1);
  94. std::cout << "Status has been set" << std::endl;
  95. std::cout << "Now I know what to read, Master!" << std::endl;
  96. if (!femtoReader->chain())
  97. {
  98. std::cout << "No chain has been found." << std::endl;
  99. }
  100. Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  101. std::cout << "eventsInTree: " << eventsInTree << std::endl;
  102. Long64_t events2read = femtoReader->chain()->GetEntries();
  103. std::cout << "Number of events to read: " << events2read << std::endl;
  104. //Read ReCentering <Q> vectors
  105. TFile *fiReCent = new TFile(reCentFile, "read");
  106. std::map<Double_t, Double_t> p_Qx_FULL_EP[fNharmonics][fArraySize][fNCent];
  107. std::map<Double_t, Double_t> p_Qy_FULL_EP[fNharmonics][fArraySize][fNCent];
  108. std::map<Double_t, Double_t> p_Qx_FULL_SP[fNharmonics][fArraySize][fNCent];
  109. std::map<Double_t, Double_t> p_Qy_FULL_SP[fNharmonics][fArraySize][fNCent];
  110. std::map<Double_t, Double_t> p_Qx_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  111. std::map<Double_t, Double_t> p_Qy_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  112. std::map<Double_t, Double_t> p_Qx_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  113. std::map<Double_t, Double_t> p_Qy_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  114. std::map<Double_t, Double_t> p_Qx_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  115. std::map<Double_t, Double_t> p_Qy_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  116. std::map<Double_t, Double_t> p_Qx_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  117. std::map<Double_t, Double_t> p_Qy_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  118. TProfile2D *prof;
  119. double binCont;
  120. std::cout << "Reading file for TPC recentering..." << std::endl;
  121. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  122. {
  123. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  124. {
  125. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
  126. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%ix_FULL_EP%i", iHarm, iVtx));
  127. for (int i = 0; i < fNBins; i++)
  128. {
  129. for (int iCent = 0; iCent < fNCent; iCent++)
  130. {
  131. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  132. if (binCont == 0.)
  133. continue;
  134. p_Qx_FULL_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  135. }
  136. }
  137. delete prof;
  138. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%iy_FULL_EP%i", iHarm, iVtx));
  139. for (int i = 0; i < fNBins; i++)
  140. {
  141. for (int iCent = 0; iCent < fNCent; iCent++)
  142. {
  143. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  144. if (binCont == 0.)
  145. continue;
  146. p_Qy_FULL_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  147. }
  148. }
  149. delete prof;
  150. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%ix_FULL_SP%i", iHarm, iVtx));
  151. for (int i = 0; i < fNBins; i++)
  152. {
  153. for (int iCent = 0; iCent < fNCent; iCent++)
  154. {
  155. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  156. if (binCont == 0.)
  157. continue;
  158. p_Qx_FULL_SP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  159. }
  160. }
  161. delete prof;
  162. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%iy_FULL_SP%i", iHarm, iVtx));
  163. for (int i = 0; i < fNBins; i++)
  164. {
  165. for (int iCent = 0; iCent < fNCent; iCent++)
  166. {
  167. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  168. if (binCont == 0.)
  169. continue;
  170. p_Qy_FULL_SP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  171. }
  172. }
  173. delete prof;
  174. for (int iEtaGap = 0; iEtaGap < NEtaGaps; iEtaGap++)
  175. {
  176. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%ix_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  177. for (int i = 0; i < fNBins; i++)
  178. {
  179. for (int iCent = 0; iCent < fNCent; iCent++)
  180. {
  181. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  182. if (binCont == 0.)
  183. continue;
  184. p_Qx_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  185. }
  186. }
  187. delete prof;
  188. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%iy_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  189. for (int i = 0; i < fNBins; i++)
  190. {
  191. for (int iCent = 0; iCent < fNCent; iCent++)
  192. {
  193. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  194. if (binCont == 0.)
  195. continue;
  196. p_Qy_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  197. }
  198. }
  199. delete prof;
  200. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%ix_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  201. for (int i = 0; i < fNBins; i++)
  202. {
  203. for (int iCent = 0; iCent < fNCent; iCent++)
  204. {
  205. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  206. if (binCont == 0.)
  207. continue;
  208. p_Qx_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  209. }
  210. }
  211. delete prof;
  212. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%iy_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  213. for (int i = 0; i < fNBins; i++)
  214. {
  215. for (int iCent = 0; iCent < fNCent; iCent++)
  216. {
  217. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  218. if (binCont == 0.)
  219. continue;
  220. p_Qy_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  221. }
  222. }
  223. delete prof;
  224. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%ix_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  225. for (int i = 0; i < fNBins; i++)
  226. {
  227. for (int iCent = 0; iCent < fNCent; iCent++)
  228. {
  229. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  230. if (binCont == 0.)
  231. continue;
  232. p_Qx_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  233. }
  234. }
  235. delete prof;
  236. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%iy_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  237. for (int i = 0; i < fNBins; i++)
  238. {
  239. for (int iCent = 0; iCent < fNCent; iCent++)
  240. {
  241. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  242. if (binCont == 0.)
  243. continue;
  244. p_Qy_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  245. }
  246. }
  247. delete prof;
  248. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%ix_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  249. for (int i = 0; i < fNBins; i++)
  250. {
  251. for (int iCent = 0; iCent < fNCent; iCent++)
  252. {
  253. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  254. if (binCont == 0.)
  255. continue;
  256. p_Qx_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  257. }
  258. }
  259. delete prof;
  260. prof = (TProfile2D *) fiReCent->Get(Form("p_Q%iy_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  261. for (int i = 0; i < fNBins; i++)
  262. {
  263. for (int iCent = 0; iCent < fNCent; iCent++)
  264. {
  265. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  266. if (binCont == 0.)
  267. continue;
  268. p_Qy_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t) (runId_min + i)] = binCont;
  269. }
  270. }
  271. delete prof;
  272. }
  273. }
  274. }
  275. //Read Shift coefficients
  276. TFile *fiShift = new TFile(shiftFile, "read");
  277. std::map<Double_t, Double_t> p_Cos2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  278. std::map<Double_t, Double_t> p_Sin2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  279. std::map<Double_t, Double_t> p_Cos2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  280. std::map<Double_t, Double_t> p_Sin2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  281. std::map<Double_t, Double_t> p_Cos2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  282. std::map<Double_t, Double_t> p_Sin2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  283. std::map<Double_t, Double_t> p_Cos2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  284. std::map<Double_t, Double_t> p_Sin2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  285. std::map<Double_t, Double_t> p_Cos2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  286. std::map<Double_t, Double_t> p_Sin2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  287. std::map<Double_t, Double_t> p_Cos2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  288. std::map<Double_t, Double_t> p_Sin2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  289. std::map<Double_t, Double_t> p_Cos3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  290. std::map<Double_t, Double_t> p_Sin3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  291. std::map<Double_t, Double_t> p_Cos3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  292. std::map<Double_t, Double_t> p_Sin3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  293. std::map<Double_t, Double_t> p_Cos3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  294. std::map<Double_t, Double_t> p_Sin3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  295. std::map<Double_t, Double_t> p_Cos3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  296. std::map<Double_t, Double_t> p_Sin3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  297. std::map<Double_t, Double_t> p_Cos3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  298. std::map<Double_t, Double_t> p_Sin3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  299. std::map<Double_t, Double_t> p_Cos3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  300. std::map<Double_t, Double_t> p_Sin3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  301. std::cout << "Reading file for TPC flattening..." << std::endl;
  302. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  303. {
  304. for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
  305. {
  306. std::cout << "iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
  307. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
  308. for (int i = 0; i < fNBins; i++)
  309. {
  310. for (int iCent = 0; iCent < fNCent; iCent++)
  311. {
  312. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  313. if (binCont == 0.)
  314. continue;
  315. p_Cos2_FULL_EP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  316. }
  317. }
  318. delete prof;
  319. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
  320. for (int i = 0; i < fNBins; i++)
  321. {
  322. for (int iCent = 0; iCent < fNCent; iCent++)
  323. {
  324. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  325. if (binCont == 0.)
  326. continue;
  327. p_Sin2_FULL_EP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  328. }
  329. }
  330. delete prof;
  331. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
  332. for (int i = 0; i < fNBins; i++)
  333. {
  334. for (int iCent = 0; iCent < fNCent; iCent++)
  335. {
  336. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  337. if (binCont == 0.)
  338. continue;
  339. p_Cos2_FULL_SP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  340. }
  341. }
  342. delete prof;
  343. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
  344. for (int i = 0; i < fNBins; i++)
  345. {
  346. for (int iCent = 0; iCent < fNCent; iCent++)
  347. {
  348. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  349. if (binCont == 0.)
  350. continue;
  351. p_Sin2_FULL_SP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  352. }
  353. }
  354. delete prof;
  355. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
  356. for (int i = 0; i < fNBins; i++)
  357. {
  358. for (int iCent = 0; iCent < fNCent; iCent++)
  359. {
  360. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  361. if (binCont == 0.)
  362. continue;
  363. p_Cos3_FULL_EP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  364. }
  365. }
  366. delete prof;
  367. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
  368. for (int i = 0; i < fNBins; i++)
  369. {
  370. for (int iCent = 0; iCent < fNCent; iCent++)
  371. {
  372. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  373. if (binCont == 0.)
  374. continue;
  375. p_Sin3_FULL_EP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  376. }
  377. }
  378. delete prof;
  379. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
  380. for (int i = 0; i < fNBins; i++)
  381. {
  382. for (int iCent = 0; iCent < fNCent; iCent++)
  383. {
  384. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  385. if (binCont == 0.)
  386. continue;
  387. p_Cos3_FULL_SP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  388. }
  389. }
  390. delete prof;
  391. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
  392. for (int i = 0; i < fNBins; i++)
  393. {
  394. for (int iCent = 0; iCent < fNCent; iCent++)
  395. {
  396. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  397. if (binCont == 0.)
  398. continue;
  399. p_Sin3_FULL_SP[iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  400. }
  401. }
  402. delete prof;
  403. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  404. {
  405. // for 2nd harmonics
  406. // EAST
  407. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  408. for (int i = 0; i < fNBins; i++)
  409. {
  410. for (int iCent = 0; iCent < fNCent; iCent++)
  411. {
  412. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  413. if (binCont == 0.)
  414. continue;
  415. p_Cos2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  416. }
  417. }
  418. delete prof;
  419. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  420. for (int i = 0; i < fNBins; i++)
  421. {
  422. for (int iCent = 0; iCent < fNCent; iCent++)
  423. {
  424. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  425. if (binCont == 0.)
  426. continue;
  427. p_Sin2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  428. }
  429. }
  430. delete prof;
  431. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  432. for (int i = 0; i < fNBins; i++)
  433. {
  434. for (int iCent = 0; iCent < fNCent; iCent++)
  435. {
  436. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  437. if (binCont == 0.)
  438. continue;
  439. p_Cos2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  440. }
  441. }
  442. delete prof;
  443. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  444. for (int i = 0; i < fNBins; i++)
  445. {
  446. for (int iCent = 0; iCent < fNCent; iCent++)
  447. {
  448. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  449. if (binCont == 0.)
  450. continue;
  451. p_Sin2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  452. }
  453. }
  454. delete prof;
  455. // WEST
  456. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  457. for (int i = 0; i < fNBins; i++)
  458. {
  459. for (int iCent = 0; iCent < fNCent; iCent++)
  460. {
  461. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  462. if (binCont == 0.)
  463. continue;
  464. p_Cos2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  465. }
  466. }
  467. delete prof;
  468. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  469. for (int i = 0; i < fNBins; i++)
  470. {
  471. for (int iCent = 0; iCent < fNCent; iCent++)
  472. {
  473. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  474. if (binCont == 0.)
  475. continue;
  476. p_Sin2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  477. }
  478. }
  479. delete prof;
  480. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  481. for (int i = 0; i < fNBins; i++)
  482. {
  483. for (int iCent = 0; iCent < fNCent; iCent++)
  484. {
  485. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  486. if (binCont == 0.)
  487. continue;
  488. p_Cos2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  489. }
  490. }
  491. delete prof;
  492. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  493. for (int i = 0; i < fNBins; i++)
  494. {
  495. for (int iCent = 0; iCent < fNCent; iCent++)
  496. {
  497. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  498. if (binCont == 0.)
  499. continue;
  500. p_Sin2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  501. }
  502. }
  503. delete prof;
  504. // for 3rd harmonics
  505. // EAST
  506. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  507. for (int i = 0; i < fNBins; i++)
  508. {
  509. for (int iCent = 0; iCent < fNCent; iCent++)
  510. {
  511. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  512. if (binCont == 0.)
  513. continue;
  514. p_Cos3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  515. }
  516. }
  517. delete prof;
  518. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  519. for (int i = 0; i < fNBins; i++)
  520. {
  521. for (int iCent = 0; iCent < fNCent; iCent++)
  522. {
  523. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  524. if (binCont == 0.)
  525. continue;
  526. p_Sin3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  527. }
  528. }
  529. delete prof;
  530. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  531. for (int i = 0; i < fNBins; i++)
  532. {
  533. for (int iCent = 0; iCent < fNCent; iCent++)
  534. {
  535. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  536. if (binCont == 0.)
  537. continue;
  538. p_Cos3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  539. }
  540. }
  541. delete prof;
  542. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  543. for (int i = 0; i < fNBins; i++)
  544. {
  545. for (int iCent = 0; iCent < fNCent; iCent++)
  546. {
  547. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  548. if (binCont == 0.)
  549. continue;
  550. p_Sin3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  551. }
  552. }
  553. delete prof;
  554. // WEST
  555. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  556. for (int i = 0; i < fNBins; i++)
  557. {
  558. for (int iCent = 0; iCent < fNCent; iCent++)
  559. {
  560. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  561. if (binCont == 0.)
  562. continue;
  563. p_Cos3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  564. }
  565. }
  566. delete prof;
  567. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  568. for (int i = 0; i < fNBins; i++)
  569. {
  570. for (int iCent = 0; iCent < fNCent; iCent++)
  571. {
  572. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  573. if (binCont == 0.)
  574. continue;
  575. p_Sin3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  576. }
  577. }
  578. delete prof;
  579. prof = (TProfile2D *) fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  580. for (int i = 0; i < fNBins; i++)
  581. {
  582. for (int iCent = 0; iCent < fNCent; iCent++)
  583. {
  584. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  585. if (binCont == 0.)
  586. continue;
  587. p_Cos3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  588. }
  589. }
  590. delete prof;
  591. prof = (TProfile2D *) fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  592. for (int i = 0; i < fNBins; i++)
  593. {
  594. for (int iCent = 0; iCent < fNCent; iCent++)
  595. {
  596. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  597. if (binCont == 0.)
  598. continue;
  599. p_Sin3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t) (runId_min + i)] = binCont;
  600. }
  601. }
  602. delete prof;
  603. }
  604. }
  605. }
  606. //Read resolution
  607. TFile *fiRes = new TFile(resFile, "read");
  608. TProfile *p_res2_EP[NEtaGaps];
  609. TProfile *p_res3_EP[NEtaGaps];
  610. TProfile *p_res2_SP[NEtaGaps];
  611. TProfile *p_res3_SP[NEtaGaps];
  612. // TProfile *p_resBBC_EP[fNharmonics+1];
  613. TProfile *p_res2_EP_3sub_AB;
  614. TProfile *p_res2_EP_3sub_AC;
  615. TProfile *p_res2_EP_3sub_BC;
  616. p_res2_EP_3sub_AB = (TProfile *) fiRes->Get("p_res2_EP_3sub_AB");
  617. p_res2_EP_3sub_AC = (TProfile *) fiRes->Get("p_res2_EP_3sub_AC");
  618. p_res2_EP_3sub_BC = (TProfile *) fiRes->Get("p_res2_EP_3sub_BC");
  619. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  620. {
  621. p_res2_EP[iGap] = (TProfile *) fiRes->Get(Form("p_res2_EP_gap%i", iGap));
  622. p_res3_EP[iGap] = (TProfile *) fiRes->Get(Form("p_res3_EP_gap%i", iGap));
  623. p_res2_SP[iGap] = (TProfile *) fiRes->Get(Form("p_res2_SP_gap%i", iGap));
  624. p_res3_SP[iGap] = (TProfile *) fiRes->Get(Form("p_res3_SP_gap%i", iGap));
  625. }
  626. // for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  627. // {
  628. // p_resBBC_EP[iHarm] = (TProfile*) fiRes->Get(Form("p_resBBC_EP%i", iHarm));
  629. // }
  630. //Read resolution fit
  631. TFile *fiResFit = new TFile(resFitFile, "read");
  632. TH1D *hResBBCFULL[fNharmonics + 1];
  633. TH1D *hResBBCHalf[fNharmonics + 1];
  634. TH1D *hResZDCFULL[fNharmonics + 1];
  635. TH1D *hResZDCHalf[fNharmonics + 1];
  636. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  637. {
  638. hResBBCFULL[iHarm] = (TH1D *) fiResFit->Get(Form("hResBBCFULL%i", iHarm));
  639. hResBBCHalf[iHarm] = (TH1D *) fiResFit->Get(Form("hResBBCHalf%i", iHarm));
  640. }
  641. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  642. {
  643. hResZDCFULL[iHarm] = (TH1D *) fiResFit->Get(Form("hResZDCFULL%i", iHarm));
  644. hResZDCHalf[iHarm] = (TH1D *) fiResFit->Get(Form("hResZDCHalf%i", iHarm));
  645. }
  646. //Manage gain correction
  647. TProfile2D *BBCgain_East[fArraySize];
  648. TProfile2D *BBCgain_West[fArraySize];
  649. TFile *fiGain = new TFile(gainBBC, "read");
  650. fiGain->cd();
  651. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  652. {
  653. BBCgain_East[iVtx] = (TProfile2D *) fiGain->Get(Form("p_gainBBC_East_Vtx%i", iVtx));
  654. BBCgain_West[iVtx] = (TProfile2D *) fiGain->Get(Form("p_gainBBC_West_Vtx%i", iVtx));
  655. }
  656. //Manage recentering BBC
  657. std::map<Double_t, Double_t> p_BBC_Qx_Full_EP[fNharmonics][fArraySize][fNCent];
  658. std::map<Double_t, Double_t> p_BBC_Qy_Full_EP[fNharmonics][fArraySize][fNCent];
  659. std::map<Double_t, Double_t> p_BBC_Qx_East_EP[fNharmonics][fArraySize][fNCent];
  660. std::map<Double_t, Double_t> p_BBC_Qy_East_EP[fNharmonics][fArraySize][fNCent];
  661. std::map<Double_t, Double_t> p_BBC_Qx_West_EP[fNharmonics][fArraySize][fNCent];
  662. std::map<Double_t, Double_t> p_BBC_Qy_West_EP[fNharmonics][fArraySize][fNCent];
  663. TFile *fiReCentBBC = new TFile(reCentFileBBC, "read");
  664. std::cout << "Reading file for BBC recentering..." << std::endl;
  665. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  666. {
  667. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  668. {
  669. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
  670. prof = (TProfile2D *) fiReCentBBC->Get(Form("p_BBC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
  671. for (int i = 0; i < fNBins; i++)
  672. {
  673. for (int iCent = 0; iCent < fNCent; iCent++)
  674. {
  675. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  676. if (binCont == 0.)
  677. continue;
  678. p_BBC_Qx_Full_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  679. }
  680. }
  681. delete prof;
  682. prof = (TProfile2D *) fiReCentBBC->Get(Form("p_BBC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
  683. for (int i = 0; i < fNBins; i++)
  684. {
  685. for (int iCent = 0; iCent < fNCent; iCent++)
  686. {
  687. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  688. if (binCont == 0.)
  689. continue;
  690. p_BBC_Qy_Full_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  691. }
  692. }
  693. delete prof;
  694. prof = (TProfile2D *) fiReCentBBC->Get(Form("p_BBC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
  695. for (int i = 0; i < fNBins; i++)
  696. {
  697. for (int iCent = 0; iCent < fNCent; iCent++)
  698. {
  699. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  700. if (binCont == 0.)
  701. continue;
  702. p_BBC_Qx_East_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  703. }
  704. }
  705. delete prof;
  706. prof = (TProfile2D *) fiReCentBBC->Get(Form("p_BBC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
  707. for (int i = 0; i < fNBins; i++)
  708. {
  709. for (int iCent = 0; iCent < fNCent; iCent++)
  710. {
  711. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  712. if (binCont == 0.)
  713. continue;
  714. p_BBC_Qy_East_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  715. }
  716. }
  717. delete prof;
  718. prof = (TProfile2D *) fiReCentBBC->Get(Form("p_BBC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
  719. for (int i = 0; i < fNBins; i++)
  720. {
  721. for (int iCent = 0; iCent < fNCent; iCent++)
  722. {
  723. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  724. if (binCont == 0.)
  725. continue;
  726. p_BBC_Qx_West_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  727. }
  728. }
  729. delete prof;
  730. prof = (TProfile2D *) fiReCentBBC->Get(Form("p_BBC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
  731. for (int i = 0; i < fNBins; i++)
  732. {
  733. for (int iCent = 0; iCent < fNCent; iCent++)
  734. {
  735. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  736. if (binCont == 0.)
  737. continue;
  738. p_BBC_Qy_West_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  739. }
  740. }
  741. delete prof;
  742. }
  743. }
  744. //Manage ShiftCorr BBC
  745. std::map<Double_t, Double_t> p_BBC_Cos_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  746. std::map<Double_t, Double_t> p_BBC_Sin_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  747. std::map<Double_t, Double_t> p_BBC_Cos_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  748. std::map<Double_t, Double_t> p_BBC_Sin_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  749. std::map<Double_t, Double_t> p_BBC_Cos_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  750. std::map<Double_t, Double_t> p_BBC_Sin_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  751. TFile *fiShiftBBC = new TFile(shiftFileBBC, "read");
  752. fiShiftBBC->cd();
  753. std::cout << "Reading file for BBC shift corrections..." << std::endl;
  754. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  755. {
  756. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  757. {
  758. for (int iShift = 0; iShift < NShiftOrderMax * 3; iShift++)
  759. {
  760. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
  761. prof = (TProfile2D *) fiShiftBBC->Get(Form("p_BBC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
  762. for (int i = 0; i < fNBins; i++)
  763. {
  764. for (int iCent = 0; iCent < fNCent; iCent++)
  765. {
  766. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  767. if (binCont == 0.)
  768. continue;
  769. p_BBC_Cos_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  770. }
  771. }
  772. delete prof;
  773. prof = (TProfile2D *) fiShiftBBC->Get(Form("p_BBC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
  774. for (int i = 0; i < fNBins; i++)
  775. {
  776. for (int iCent = 0; iCent < fNCent; iCent++)
  777. {
  778. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  779. if (binCont == 0.)
  780. continue;
  781. p_BBC_Sin_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  782. }
  783. }
  784. delete prof;
  785. prof = (TProfile2D *) fiShiftBBC->Get(Form("p_BBC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
  786. for (int i = 0; i < fNBins; i++)
  787. {
  788. for (int iCent = 0; iCent < fNCent; iCent++)
  789. {
  790. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  791. if (binCont == 0.)
  792. continue;
  793. p_BBC_Cos_East_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  794. }
  795. }
  796. delete prof;
  797. prof = (TProfile2D *) fiShiftBBC->Get(Form("p_BBC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
  798. for (int i = 0; i < fNBins; i++)
  799. {
  800. for (int iCent = 0; iCent < fNCent; iCent++)
  801. {
  802. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  803. if (binCont == 0.)
  804. continue;
  805. p_BBC_Sin_East_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  806. }
  807. }
  808. delete prof;
  809. prof = (TProfile2D *) fiShiftBBC->Get(Form("p_BBC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
  810. for (int i = 0; i < fNBins; i++)
  811. {
  812. for (int iCent = 0; iCent < fNCent; iCent++)
  813. {
  814. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  815. if (binCont == 0.)
  816. continue;
  817. p_BBC_Cos_West_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  818. }
  819. }
  820. delete prof;
  821. prof = (TProfile2D *) fiShiftBBC->Get(Form("p_BBC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
  822. for (int i = 0; i < fNBins; i++)
  823. {
  824. for (int iCent = 0; iCent < fNCent; iCent++)
  825. {
  826. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  827. if (binCont == 0.)
  828. continue;
  829. p_BBC_Sin_West_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  830. }
  831. }
  832. delete prof;
  833. }
  834. }
  835. }
  836. //Manage recentering ZDC
  837. std::map<Double_t, Double_t> p_ZDC_Qx_Full_EP[fNharmonics][fArraySize][fNCent];
  838. std::map<Double_t, Double_t> p_ZDC_Qy_Full_EP[fNharmonics][fArraySize][fNCent];
  839. std::map<Double_t, Double_t> p_ZDC_Qx_East_EP[fNharmonics][fArraySize][fNCent];
  840. std::map<Double_t, Double_t> p_ZDC_Qy_East_EP[fNharmonics][fArraySize][fNCent];
  841. std::map<Double_t, Double_t> p_ZDC_Qx_West_EP[fNharmonics][fArraySize][fNCent];
  842. std::map<Double_t, Double_t> p_ZDC_Qy_West_EP[fNharmonics][fArraySize][fNCent];
  843. TFile *fiReCentZDC = new TFile(reCentFileZDC, "read");
  844. std::cout << "Reading file for ZDC recentering..." << std::endl;
  845. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  846. {
  847. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  848. {
  849. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
  850. prof = (TProfile2D *) fiReCentZDC->Get(Form("p_ZDC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
  851. for (int i = 0; i < fNBins; i++)
  852. {
  853. for (int iCent = 0; iCent < fNCent; iCent++)
  854. {
  855. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  856. if (binCont == 0.)
  857. continue;
  858. p_ZDC_Qx_Full_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  859. }
  860. }
  861. delete prof;
  862. prof = (TProfile2D *) fiReCentZDC->Get(Form("p_ZDC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
  863. for (int i = 0; i < fNBins; i++)
  864. {
  865. for (int iCent = 0; iCent < fNCent; iCent++)
  866. {
  867. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  868. if (binCont == 0.)
  869. continue;
  870. p_ZDC_Qy_Full_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  871. }
  872. }
  873. delete prof;
  874. prof = (TProfile2D *) fiReCentZDC->Get(Form("p_ZDC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
  875. for (int i = 0; i < fNBins; i++)
  876. {
  877. for (int iCent = 0; iCent < fNCent; iCent++)
  878. {
  879. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  880. if (binCont == 0.)
  881. continue;
  882. p_ZDC_Qx_East_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  883. }
  884. }
  885. delete prof;
  886. prof = (TProfile2D *) fiReCentZDC->Get(Form("p_ZDC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
  887. for (int i = 0; i < fNBins; i++)
  888. {
  889. for (int iCent = 0; iCent < fNCent; iCent++)
  890. {
  891. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  892. if (binCont == 0.)
  893. continue;
  894. p_ZDC_Qy_East_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  895. }
  896. }
  897. delete prof;
  898. prof = (TProfile2D *) fiReCentZDC->Get(Form("p_ZDC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
  899. for (int i = 0; i < fNBins; i++)
  900. {
  901. for (int iCent = 0; iCent < fNCent; iCent++)
  902. {
  903. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  904. if (binCont == 0.)
  905. continue;
  906. p_ZDC_Qx_West_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  907. }
  908. }
  909. delete prof;
  910. prof = (TProfile2D *) fiReCentZDC->Get(Form("p_ZDC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
  911. for (int i = 0; i < fNBins; i++)
  912. {
  913. for (int iCent = 0; iCent < fNCent; iCent++)
  914. {
  915. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  916. if (binCont == 0.)
  917. continue;
  918. p_ZDC_Qy_West_EP[iHarm][iVtx][iCent][(Double_t) (runId_min + i)] = binCont;
  919. }
  920. }
  921. delete prof;
  922. }
  923. }
  924. //Manage ShiftCorr ZDC
  925. std::map<Double_t, Double_t> p_ZDC_Cos_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  926. std::map<Double_t, Double_t> p_ZDC_Sin_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  927. std::map<Double_t, Double_t> p_ZDC_Cos_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  928. std::map<Double_t, Double_t> p_ZDC_Sin_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  929. std::map<Double_t, Double_t> p_ZDC_Cos_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  930. std::map<Double_t, Double_t> p_ZDC_Sin_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  931. TFile *fiShiftZDC = new TFile(shiftFileZDC, "read");
  932. fiShiftZDC->cd();
  933. std::cout << "Reading file for ZDC shift corrections..." << std::endl;
  934. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  935. {
  936. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  937. {
  938. for (int iShift = 0; iShift < NShiftOrderMax * 3; iShift++)
  939. {
  940. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
  941. prof = (TProfile2D *) fiShiftZDC->Get(Form("p_ZDC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
  942. for (int i = 0; i < fNBins; i++)
  943. {
  944. for (int iCent = 0; iCent < fNCent; iCent++)
  945. {
  946. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  947. if (binCont == 0.)
  948. continue;
  949. p_ZDC_Cos_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  950. }
  951. }
  952. delete prof;
  953. prof = (TProfile2D *) fiShiftZDC->Get(Form("p_ZDC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
  954. for (int i = 0; i < fNBins; i++)
  955. {
  956. for (int iCent = 0; iCent < fNCent; iCent++)
  957. {
  958. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  959. if (binCont == 0.)
  960. continue;
  961. p_ZDC_Sin_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  962. }
  963. }
  964. delete prof;
  965. prof = (TProfile2D *) fiShiftZDC->Get(Form("p_ZDC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
  966. for (int i = 0; i < fNBins; i++)
  967. {
  968. for (int iCent = 0; iCent < fNCent; iCent++)
  969. {
  970. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  971. if (binCont == 0.)
  972. continue;
  973. p_ZDC_Cos_East_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  974. }
  975. }
  976. delete prof;
  977. prof = (TProfile2D *) fiShiftZDC->Get(Form("p_ZDC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
  978. for (int i = 0; i < fNBins; i++)
  979. {
  980. for (int iCent = 0; iCent < fNCent; iCent++)
  981. {
  982. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  983. if (binCont == 0.)
  984. continue;
  985. p_ZDC_Sin_East_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  986. }
  987. }
  988. delete prof;
  989. prof = (TProfile2D *) fiShiftZDC->Get(Form("p_ZDC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
  990. for (int i = 0; i < fNBins; i++)
  991. {
  992. for (int iCent = 0; iCent < fNCent; iCent++)
  993. {
  994. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  995. if (binCont == 0.)
  996. continue;
  997. p_ZDC_Cos_West_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  998. }
  999. }
  1000. delete prof;
  1001. prof = (TProfile2D *) fiShiftZDC->Get(Form("p_ZDC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
  1002. for (int i = 0; i < fNBins; i++)
  1003. {
  1004. for (int iCent = 0; iCent < fNCent; iCent++)
  1005. {
  1006. binCont = prof->GetBinContent(prof->FindBin((Double_t) (runId_min + i), (Double_t) iCent));
  1007. if (binCont == 0.)
  1008. continue;
  1009. p_ZDC_Sin_West_EP[iHarm][iVtx][iShift][iCent][(Double_t) (runId_min + i)] = binCont;
  1010. }
  1011. }
  1012. delete prof;
  1013. }
  1014. }
  1015. }
  1016. //Manage output TProfiles
  1017. TProfile2D *p_v2run_EP[3][NEtaGaps][16]; // runId vs pt
  1018. TProfile2D *p_v3run_EP[3][NEtaGaps][16]; // runId vs pt
  1019. TProfile2D *p_v2run_BBC[3][16];
  1020. TProfile2D *p_v2run_ZDC[3][16];
  1021. TProfile2D *p_v2run_ZDC_East[3][16];
  1022. TProfile2D *p_v2run_ZDC_West[3][16];
  1023. // TProfile2D *p_v2_sys_Nhits_EP[3][NEtaGaps][16][2]; // runId vs pt
  1024. // TProfile2D *p_v3_sys_Nhits_EP[3][NEtaGaps][16][2]; // runId vs pt
  1025. // TProfile2D *p_v2_sys_DCA_EP[3][NEtaGaps][16][2]; // runId vs pt
  1026. // TProfile2D *p_v3_sys_DCA_EP[3][NEtaGaps][16][2]; // runId vs pt
  1027. // TProfile2D *p_v2_sys_Eta_EP[3][NEtaGaps][16]; // runId vs pt
  1028. // TProfile2D *p_v3_sys_Eta_EP[3][NEtaGaps][16]; // runId vs pt
  1029. for (int iPID = 0; iPID < 3; iPID++)
  1030. {
  1031. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  1032. {
  1033. for (int iCent = 0; iCent < 16; iCent++)
  1034. {
  1035. p_v2run_EP[iPID][iGap][iCent] = new TProfile2D(Form("p_v2_run_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1036. Form("p_v2_run_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1037. fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1038. p_v3run_EP[iPID][iGap][iCent] = new TProfile2D(Form("p_v3_run_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1039. Form("p_v3_run_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1040. fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1041. // p_v2_sys_Eta_EP[iPID][iGap][iCent] = new TProfile2D(Form("p_v2_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1042. // Form("p_v2_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1043. // fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1044. // p_v3_sys_Eta_EP[iPID][iGap][iCent] = new TProfile2D(Form("p_v3_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1045. // Form("p_v3_sys_Eta_EP_particle%i_gap%i_cent%i", iPID, iGap, iCent),
  1046. // fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1047. // for (int iCut = 0; iCut < 2; iCut++)
  1048. // {
  1049. // p_v2_sys_Nhits_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
  1050. // Form("p_v2_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
  1051. // Form("p_v2_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
  1052. // runId_max + 0.5, 100, 0.15, 5.15);
  1053. // p_v3_sys_Nhits_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
  1054. // Form("p_v3_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
  1055. // Form("p_v3_sys_Nhits_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
  1056. // runId_max + 0.5, 100, 0.15, 5.15);
  1057. //
  1058. // p_v2_sys_DCA_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
  1059. // Form("p_v2_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
  1060. // Form("p_v2_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
  1061. // runId_max + 0.5, 100, 0.15, 5.15);
  1062. // p_v3_sys_DCA_EP[iPID][iGap][iCent][iCut] = new TProfile2D(
  1063. // Form("p_v3_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut),
  1064. // Form("p_v3_sys_DCA_EP_particle%i_gap%i_cent%i_sys%i", iPID, iGap, iCent, iCut), fNBins, runId_min - 0.5,
  1065. // runId_max + 0.5, 100, 0.15, 5.15);
  1066. // }
  1067. if (iGap == 0)
  1068. {
  1069. p_v2run_BBC[iPID][iCent] = new TProfile2D(Form("p_v2_run_BBC_particle%i_cent%i", iPID, iCent),
  1070. Form("p_v2_run_BBC_particle%i_cent%i", iPID, iCent), fNBins,
  1071. runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1072. p_v2run_ZDC[iPID][iCent] = new TProfile2D(Form("p_v2_run_ZDC_particle%i_cent%i", iPID, iCent),
  1073. Form("p_v2_run_ZDC_particle%i_cent%i", iPID, iCent), fNBins,
  1074. runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1075. p_v2run_ZDC_East[iPID][iCent] = new TProfile2D(Form("p_v2_run_ZDC_East_particle%i_cent%i", iPID, iCent),
  1076. Form("p_v2_run_ZDC_East_particle%i_cent%i", iPID, iCent),
  1077. fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1078. p_v2run_ZDC_West[iPID][iCent] = new TProfile2D(Form("p_v2_run_ZDC_West_particle%i_cent%i", iPID, iCent),
  1079. Form("p_v2_run_ZDC_West_particle%i_cent%i", iPID, iCent),
  1080. fNBins, runId_min - 0.5, runId_max + 0.5, 100, 0.15, 5.15);
  1081. }
  1082. }
  1083. }
  1084. }
  1085. //Counters
  1086. Int_t fUsedTracks = 0;
  1087. Int_t fQcountFull = 0;
  1088. Int_t fQcountFullEast = 0;
  1089. Int_t fQcountFullWest = 0;
  1090. Int_t fQcountWest[NEtaGaps];
  1091. Int_t fQcountEast[NEtaGaps];
  1092. //Q-vectors
  1093. TVector2 fQv2FullEP;
  1094. TVector2 fQv2FullSP;
  1095. TVector2 fQv3FullEP;
  1096. TVector2 fQv3FullSP;
  1097. TVector2 fQv2WestEP[NEtaGaps];
  1098. TVector2 fQv2EastEP[NEtaGaps];
  1099. TVector2 fQv2WestSP[NEtaGaps];
  1100. TVector2 fQv2EastSP[NEtaGaps];
  1101. TVector2 fQv3WestEP[NEtaGaps];
  1102. TVector2 fQv3EastEP[NEtaGaps];
  1103. TVector2 fQv3WestSP[NEtaGaps];
  1104. TVector2 fQv3EastSP[NEtaGaps];
  1105. TVector2 QvTMP;
  1106. Double_t weight;
  1107. Int_t VtxSign;
  1108. Double_t psiShiftedTMP;
  1109. Double_t PsiSin;
  1110. Double_t PsiCos;
  1111. Long64_t goodEventCounter = 0;
  1112. Double_t Psi2FullEP, Psi3FullEP;
  1113. Double_t Psi2EastEP[NEtaGaps], Psi3EastEP[NEtaGaps];
  1114. Double_t Psi2WestEP[NEtaGaps], Psi3WestEP[NEtaGaps];
  1115. Double_t Psi_ReCenter, dPsi;
  1116. Int_t bin;
  1117. Double_t res2EP, res3EP, res2EPsqr, res3EPsqr;
  1118. Double_t flow2EP;
  1119. Double_t flow3EP;
  1120. Int_t i_part;
  1121. Double_t flow2EPBBC, res2EPBBC;
  1122. Double_t psiEP_rec, psiEP_BBC_Full[fNharmonics + 1], shiftCos, shiftSin;
  1123. Double_t psiEP_BBC_East[fNharmonics + 1], psiEP_BBC_West[fNharmonics + 1];
  1124. // ZDC related EP determination
  1125. Double_t flow2EPZDC, res2EPZDC;
  1126. Double_t psiEP_ZDC_Full[fNharmonics + 1], psiEP_ZDC_East[fNharmonics + 1], psiEP_ZDC_West[fNharmonics + 1];
  1127. TVector2 qRecEastZDC[fNharmonics + 1], qRecWestZDC[fNharmonics + 1], qRecFullZDC[fNharmonics + 1];
  1128. Double_t qxEast, qyEast, adcxEast, adcyEast;
  1129. Double_t qxWest, qyWest, adcxWest, adcyWest;
  1130. Double_t qxFull, qyFull, adcxFull, adcyFull;
  1131. Double_t qxRecEast, qyRecEast;
  1132. Double_t qxRecWest, qyRecWest;
  1133. Double_t qxRecFull, qyRecFull;
  1134. Double_t resEastWestZDC[fNharmonics + 1];
  1135. std::pair<Double_t, Double_t> XY_pid;
  1136. // Loop over events
  1137. for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
  1138. {
  1139. if ((iEvent + 1) % 1000 == 0)
  1140. {
  1141. std::cout << "Working on event #[" << (iEvent + 1)
  1142. << "/" << events2read << "]" << std::endl;
  1143. }
  1144. Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
  1145. if (!readEvent)
  1146. {
  1147. std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
  1148. break;
  1149. }
  1150. // Retrieve femtoDst
  1151. StFemtoDst *dst = femtoReader->femtoDst();
  1152. // Retrieve event information
  1153. StFemtoEvent *event = dst->event();
  1154. if (!event)
  1155. {
  1156. std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
  1157. break;
  1158. }
  1159. // Event selection
  1160. if (!isGoodEvent(event)) continue;
  1161. goodEventCounter++;
  1162. TVector3 pVtx = event->primaryVertex();
  1163. // Track analysis
  1164. Int_t nTracks = dst->numberOfTracks();
  1165. //Counters
  1166. fUsedTracks = 0;
  1167. fQcountFull = 0;
  1168. fQcountFullEast = 0;
  1169. fQcountFullWest = 0;
  1170. psiShiftedTMP = 0.;
  1171. PsiSin = 0.;
  1172. PsiCos = 0.;
  1173. fQv2FullEP.Set(0., 0.);
  1174. fQv2FullSP.Set(0., 0.);
  1175. fQv3FullEP.Set(0., 0.);
  1176. fQv3FullSP.Set(0., 0.);
  1177. QvTMP.Set(0., 0.);
  1178. for (int i = 0; i < NEtaGaps; i++)
  1179. {
  1180. fQcountWest[i] = 0;
  1181. fQcountEast[i] = 0;
  1182. fQv2WestEP[i].Set(0., 0.);
  1183. fQv2EastEP[i].Set(0., 0.);
  1184. fQv2WestSP[i].Set(0., 0.);
  1185. fQv2EastSP[i].Set(0., 0.);
  1186. fQv3WestEP[i].Set(0., 0.);
  1187. fQv3EastEP[i].Set(0., 0.);
  1188. fQv3WestSP[i].Set(0., 0.);
  1189. fQv3EastSP[i].Set(0., 0.);
  1190. }
  1191. // Vertex sign
  1192. VtxSign = GetVzBin(pVtx.Z());
  1193. if (VtxSign == -1) continue;
  1194. //BBC EP
  1195. Double_t adcNormEastInner =
  1196. BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),
  1197. BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 1, 6) / 6.;
  1198. Double_t adcNormEastOuter =
  1199. BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),
  1200. BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 7, 16) / 10.;
  1201. Double_t adcNormWestInner =
  1202. BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),
  1203. BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 1, 6) / 6.;
  1204. Double_t adcNormWestOuter =
  1205. BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),
  1206. BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()), 7, 16) / 10.;
  1207. float egain[16] = {
  1208. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 1)) /
  1209. adcNormEastInner),
  1210. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 2)) /
  1211. adcNormEastInner),
  1212. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 3)) /
  1213. adcNormEastInner),
  1214. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 4)) /
  1215. adcNormEastInner),
  1216. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 5)) /
  1217. adcNormEastInner),
  1218. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 6)) /
  1219. adcNormEastInner),
  1220. (float) ((0.5) * BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 7)) /
  1221. adcNormEastOuter),
  1222. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 8)) /
  1223. adcNormEastOuter),
  1224. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 9)) /
  1225. adcNormEastOuter),
  1226. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 10)) /
  1227. adcNormEastOuter),
  1228. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 11)) /
  1229. adcNormEastOuter),
  1230. (float) ((0.5) * BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 12)) /
  1231. adcNormEastOuter),
  1232. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 13)) /
  1233. adcNormEastOuter),
  1234. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 14)) /
  1235. adcNormEastOuter),
  1236. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 15)) /
  1237. adcNormEastOuter),
  1238. (float) (BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(), 16)) /
  1239. adcNormEastOuter)};
  1240. float wgain[16] = {
  1241. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 1)) /
  1242. adcNormWestInner),
  1243. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 2)) /
  1244. adcNormWestInner),
  1245. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 3)) /
  1246. adcNormWestInner),
  1247. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 4)) /
  1248. adcNormWestInner),
  1249. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 5)) /
  1250. adcNormWestInner),
  1251. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 6)) /
  1252. adcNormWestInner),
  1253. (float) ((0.5) * BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 7)) /
  1254. adcNormWestOuter),
  1255. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 8)) /
  1256. adcNormWestOuter),
  1257. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 9)) /
  1258. adcNormWestOuter),
  1259. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 10)) /
  1260. adcNormWestOuter),
  1261. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 11)) /
  1262. adcNormWestOuter),
  1263. (float) ((0.5) * BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 12)) /
  1264. adcNormWestOuter),
  1265. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 13)) /
  1266. adcNormWestOuter),
  1267. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 14)) /
  1268. adcNormWestOuter),
  1269. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 15)) /
  1270. adcNormWestOuter),
  1271. (float) (BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(), 16)) /
  1272. adcNormWestOuter)};
  1273. float nHitE_Gain[16], nHitW_Gain[16];
  1274. for (int i = 0; i < 16; i++)
  1275. {
  1276. nHitE_Gain[i] = event->bbcAdcEast(i) / egain[i];
  1277. nHitW_Gain[i] = event->bbcAdcWest(i) / wgain[i];
  1278. }
  1279. Double_t psiFull[fNharmonics + 1];
  1280. TVector2 qGainFull[fNharmonics + 1];
  1281. TVector2 qRecFull[fNharmonics + 1];
  1282. Double_t psiEast[fNharmonics + 1];
  1283. TVector2 qGainEast[fNharmonics + 1];
  1284. TVector2 qRecEast[fNharmonics + 1];
  1285. Double_t psiWest[fNharmonics + 1];
  1286. TVector2 qGainWest[fNharmonics + 1];
  1287. TVector2 qRecWest[fNharmonics + 1];
  1288. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  1289. {
  1290. //BBC Full
  1291. //v1*y > 0 by conventions most reasonable to the STAR event plane
  1292. //therefore the west event plane has the right sign for the event plane
  1293. //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
  1294. //The east event plane points opposite that of the west and the west finds the "correct" EP for the STAR coordinates, use West-East
  1295. Double_t QVector_GainX = 0., QVector_GainY = 0.;
  1296. Double_t EventPlaneGainFull = 0.;
  1297. Float_t eXFsum_Gain = 0., wXFsum_Gain = 0., eYFsum_Gain = 0., wYFsum_Gain = 0., eWFgt_Gain = 0., wWFgt_Gain = 0.;
  1298. TVector2 QFullGain;
  1299. for (int iTile = 0; iTile < 16; iTile++)
  1300. {
  1301. eXFsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
  1302. wXFsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];
  1303. eYFsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
  1304. wYFsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];
  1305. eWFgt_Gain += nHitE_Gain[iTile];
  1306. wWFgt_Gain += nHitW_Gain[iTile];
  1307. }
  1308. QVector_GainX = (eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wXFsum_Gain / wWFgt_Gain + eXFsum_Gain / eWFgt_Gain : 0.;
  1309. QVector_GainY = (eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wYFsum_Gain / wWFgt_Gain + eYFsum_Gain / eWFgt_Gain : 0.;
  1310. QFullGain.Set(QVector_GainX, QVector_GainY);
  1311. if (QFullGain.Mod())
  1312. {
  1313. EventPlaneGainFull = QFullGain.Phi();
  1314. if (EventPlaneGainFull < 0.0)
  1315. {
  1316. EventPlaneGainFull += 2. * TMath::Pi() / (iHarm + 1.);
  1317. }
  1318. }
  1319. psiFull[iHarm] = EventPlaneGainFull;
  1320. qGainFull[iHarm] = QFullGain;
  1321. //BBC West
  1322. //v1*y > 0 by conventions most reasonable to the STAR event plane
  1323. //therefore the west event plane has the right sign for the event plane
  1324. //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
  1325. Double_t QVectorW_GainX = 0., QVectorW_GainY = 0.;
  1326. Double_t EventPlaneGainWest = 0.;
  1327. Float_t wXsum_Gain = 0., wYsum_Gain = 0., wWgt_Gain = 0.;
  1328. TVector2 QWestGain;
  1329. for (int iTile = 0; iTile < 16; iTile++)
  1330. {
  1331. wXsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];
  1332. wYsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(1, iTile)) * nHitW_Gain[iTile];
  1333. wWgt_Gain += nHitW_Gain[iTile];
  1334. }
  1335. QVectorW_GainX = (wWgt_Gain > 0.0) ? wXsum_Gain / wWgt_Gain : 0.0;
  1336. QVectorW_GainY = (wWgt_Gain > 0.0) ? wYsum_Gain / wWgt_Gain : 0.0;
  1337. QWestGain.Set(QVectorW_GainX, QVectorW_GainY);
  1338. if (QWestGain.Mod())
  1339. {
  1340. EventPlaneGainWest = QWestGain.Phi();
  1341. if (EventPlaneGainWest < 0.0)
  1342. {
  1343. EventPlaneGainWest += 2. * TMath::Pi() / (iHarm + 1.);
  1344. }
  1345. }
  1346. psiWest[iHarm] = EventPlaneGainWest;
  1347. qGainWest[iHarm] = QWestGain;
  1348. //BBC East
  1349. //in BBC_GetPhi the fact that the tile numbering scheme is reversed about the y axis for the east BBC is accounted for
  1350. //The east event plane points opposite that of the west since the spectators are on different sides -- see west comments
  1351. Double_t QVectorE_GainX = 0., QVectorE_GainY = 0.;
  1352. Double_t EventPlaneGainEast = 0.;
  1353. Float_t eXsum_Gain = 0., eYsum_Gain = 0., eWgt_Gain = 0.;
  1354. TVector2 QEastGain;
  1355. for (int iTile = 0; iTile < 16; iTile++)
  1356. {
  1357. eXsum_Gain += cos((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
  1358. eYsum_Gain += sin((iHarm + 1.) * BBC_GetPhi(0, iTile)) * nHitE_Gain[iTile];
  1359. eWgt_Gain += nHitE_Gain[iTile];
  1360. }
  1361. QVectorE_GainX = (eWgt_Gain > 0.0) ? eXsum_Gain / eWgt_Gain : 0.0;
  1362. QVectorE_GainY = (eWgt_Gain > 0.0) ? eYsum_Gain / eWgt_Gain : 0.0;
  1363. QEastGain.Set(QVectorE_GainX, QVectorE_GainY);
  1364. if (QEastGain.Mod())
  1365. {
  1366. EventPlaneGainEast = QEastGain.Phi();
  1367. if (EventPlaneGainEast < 0.0)
  1368. {
  1369. EventPlaneGainEast += 2. * TMath::Pi() / (iHarm + 1.);
  1370. }
  1371. }
  1372. psiEast[iHarm] = EventPlaneGainEast;
  1373. qGainEast[iHarm] = QEastGain;
  1374. }
  1375. // Apply recentering
  1376. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  1377. {
  1378. TVector2 recFullTMP(
  1379. (Double_t) p_BBC_Qx_Full_EP[iHarm][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()],
  1380. (Double_t) p_BBC_Qy_Full_EP[iHarm][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1381. TVector2 recEastTMP(
  1382. (Double_t) p_BBC_Qx_East_EP[iHarm][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()],
  1383. (Double_t) p_BBC_Qy_East_EP[iHarm][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1384. TVector2 recWestTMP(
  1385. (Double_t) p_BBC_Qx_West_EP[iHarm][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()],
  1386. (Double_t) p_BBC_Qy_West_EP[iHarm][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1387. qRecEast[iHarm] = qGainEast[iHarm] - recEastTMP;
  1388. qRecWest[iHarm] = qGainWest[iHarm] - recWestTMP;
  1389. qRecFull[iHarm] = qGainFull[iHarm] - recFullTMP;
  1390. }
  1391. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  1392. {
  1393. // Full BBC
  1394. if (qRecFull[iHarm].Mod())
  1395. {
  1396. psiEP_rec = qRecFull[iHarm].Phi() / (iHarm + 1.);
  1397. if (psiEP_rec < 0.)
  1398. {
  1399. psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
  1400. }
  1401. }
  1402. dPsi = 0.;
  1403. for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
  1404. {
  1405. shiftCos = p_BBC_Cos_Full_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1406. shiftSin = p_BBC_Sin_Full_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1407. dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
  1408. (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
  1409. 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
  1410. }
  1411. psiEP_BBC_Full[iHarm] = psiEP_rec + dPsi;
  1412. if (psiEP_BBC_Full[iHarm] < 0.0)
  1413. psiEP_BBC_Full[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
  1414. if (psiEP_BBC_Full[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
  1415. psiEP_BBC_Full[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
  1416. // East BBC
  1417. if (qRecEast[iHarm].Mod())
  1418. {
  1419. psiEP_rec = qRecEast[iHarm].Phi() / (iHarm + 1.);
  1420. if (psiEP_rec < 0.)
  1421. {
  1422. psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
  1423. }
  1424. }
  1425. dPsi = 0.;
  1426. for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
  1427. {
  1428. shiftCos = p_BBC_Cos_East_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1429. shiftSin = p_BBC_Sin_East_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1430. dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
  1431. (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
  1432. 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
  1433. }
  1434. psiEP_BBC_East[iHarm] = psiEP_rec + dPsi;
  1435. if (psiEP_BBC_East[iHarm] < 0.0)
  1436. psiEP_BBC_East[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
  1437. if (psiEP_BBC_East[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
  1438. psiEP_BBC_East[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
  1439. // West BBC
  1440. if (qRecWest[iHarm].Mod())
  1441. {
  1442. psiEP_rec = qRecWest[iHarm].Phi() / (iHarm + 1.);
  1443. if (psiEP_rec < 0.)
  1444. {
  1445. psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
  1446. }
  1447. }
  1448. dPsi = 0.;
  1449. for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
  1450. {
  1451. shiftCos = p_BBC_Cos_West_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1452. shiftSin = p_BBC_Sin_West_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1453. dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
  1454. (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
  1455. 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
  1456. }
  1457. psiEP_BBC_West[iHarm] = psiEP_rec + dPsi;
  1458. if (psiEP_BBC_West[iHarm] < 0.0)
  1459. psiEP_BBC_West[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
  1460. if (psiEP_BBC_West[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
  1461. psiEP_BBC_West[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
  1462. }
  1463. //EP from BBC is ready: psiEP_BBC_Full[iHarm], where iHarm = 0 for 1st harmonics, 1 for 2nd & 2 for 3rd.
  1464. // ZDC EP
  1465. if (event->zdcSumAdcEast() < 1.)
  1466. continue;
  1467. if (event->zdcSumAdcWest() < 1.)
  1468. continue;
  1469. qxEast = 0.;
  1470. qyEast = 0.;
  1471. adcxEast = 0.;
  1472. adcyEast = 0.;
  1473. qxWest = 0.;
  1474. qyWest = 0.;
  1475. adcxWest = 0.;
  1476. adcyWest = 0.;
  1477. qxFull = 0.;
  1478. qyFull = 0.;
  1479. adcxFull = 0.;
  1480. adcyFull = 0.;
  1481. // Vertical ZDC-SMD
  1482. for (int iStrip = 0; iStrip < fNZdcSmdStripsVertical; iStrip++)
  1483. {
  1484. qxEast += event->zdcSmdEastVertical(iStrip) * GetZDCPosition(0, 0, iStrip);
  1485. adcxEast += event->zdcSmdEastVertical(iStrip);
  1486. qxWest += event->zdcSmdWestVertical(iStrip) * GetZDCPosition(1, 0, iStrip);
  1487. adcxWest += event->zdcSmdWestVertical(iStrip);
  1488. }
  1489. // Horizontal ZDC-SMD
  1490. for (int iStrip = 0; iStrip < fNZdcSmdStripsHorizontal; iStrip++)
  1491. {
  1492. qyEast += event->zdcSmdEastHorizontal(iStrip) * GetZDCPosition(0, 1, iStrip);
  1493. adcyEast += event->zdcSmdEastHorizontal(iStrip);
  1494. qyWest += event->zdcSmdWestHorizontal(iStrip) * GetZDCPosition(1, 1, iStrip);
  1495. adcyWest += event->zdcSmdWestHorizontal(iStrip);
  1496. }
  1497. if (adcxEast > 0)
  1498. qxEast /= adcxEast;
  1499. else
  1500. qxEast = -999.;
  1501. if (adcyEast > 0)
  1502. qyEast /= adcyEast;
  1503. else
  1504. qyEast = -999.;
  1505. if (adcxWest > 0)
  1506. qxWest /= adcxWest;
  1507. else
  1508. qxWest = -999.;
  1509. if (adcyWest > 0)
  1510. qyWest /= adcyWest;
  1511. else
  1512. qyWest = -999.;
  1513. // Full ZDC-SMD Q-vector
  1514. if (qxEast == -999. || qxWest == -999.)
  1515. qxFull = -999.;
  1516. else
  1517. qxFull = qxEast - qxWest;
  1518. if (qyEast == -999. || qyWest == -999.)
  1519. qyFull = -999.;
  1520. else
  1521. qyFull = qyEast - qyWest;
  1522. if (qxEast != -999.)
  1523. qxRecEast = qxEast - (Double_t) p_ZDC_Qx_East_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()];
  1524. if (qyEast != -999.)
  1525. qyRecEast = qyEast - (Double_t) p_ZDC_Qy_East_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()];
  1526. if (qxWest != -999.)
  1527. qxRecWest = qxWest - (Double_t) p_ZDC_Qx_West_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()];
  1528. if (qyWest != -999.)
  1529. qyRecWest = qyWest - (Double_t) p_ZDC_Qy_West_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()];
  1530. if (qxFull != -999.)
  1531. qxRecFull = qxFull - (Double_t) p_ZDC_Qx_Full_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()];
  1532. if (qyFull != -999.)
  1533. qyRecFull = qyFull - (Double_t) p_ZDC_Qy_Full_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()];
  1534. if (qxEast == -999.)
  1535. qxRecEast = qxEast;
  1536. if (qyEast == -999.)
  1537. qyRecEast = qyEast;
  1538. if (qxWest == -999.)
  1539. qxRecWest = qxWest;
  1540. if (qyWest == -999.)
  1541. qyRecWest = qyWest;
  1542. if (qxFull == -999.)
  1543. qxRecFull = qxFull;
  1544. if (qyFull == -999.)
  1545. qyRecFull = qyFull;
  1546. PsiSin = 0.;
  1547. PsiCos = 0.;
  1548. psiShiftedTMP = 0.;
  1549. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  1550. {
  1551. qRecEastZDC[iHarm].Set(qxRecEast, qyRecEast);
  1552. qRecWestZDC[iHarm].Set(qxRecWest, qyRecWest);
  1553. qRecFullZDC[iHarm].Set(qxRecFull, qyRecFull);
  1554. // Full ZDC
  1555. if (qRecFullZDC[iHarm].Mod())
  1556. {
  1557. psiEP_rec = qRecFullZDC[iHarm].Phi() / (iHarm + 1.);
  1558. if (psiEP_rec < 0.)
  1559. {
  1560. psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
  1561. }
  1562. }
  1563. dPsi = 0.;
  1564. for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
  1565. {
  1566. shiftCos = p_ZDC_Cos_Full_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1567. shiftSin = p_ZDC_Sin_Full_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1568. dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
  1569. (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
  1570. 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
  1571. }
  1572. psiEP_ZDC_Full[iHarm] = psiEP_rec + dPsi;
  1573. if (psiEP_ZDC_Full[iHarm] < 0.0)
  1574. psiEP_ZDC_Full[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
  1575. if (psiEP_ZDC_Full[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
  1576. psiEP_ZDC_Full[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
  1577. // East ZDC
  1578. if (qRecEastZDC[iHarm].Mod())
  1579. {
  1580. psiEP_rec = qRecEastZDC[iHarm].Phi() / (iHarm + 1.);
  1581. if (psiEP_rec < 0.)
  1582. {
  1583. psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
  1584. }
  1585. }
  1586. dPsi = 0.;
  1587. for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
  1588. {
  1589. shiftCos = p_ZDC_Cos_East_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1590. shiftSin = p_ZDC_Sin_East_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1591. dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
  1592. (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
  1593. 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
  1594. }
  1595. psiEP_ZDC_East[iHarm] = psiEP_rec + dPsi;
  1596. if (psiEP_ZDC_East[iHarm] < 0.0)
  1597. psiEP_ZDC_East[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
  1598. if (psiEP_ZDC_East[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
  1599. psiEP_ZDC_East[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
  1600. // West ZDC
  1601. if (qRecWestZDC[iHarm].Mod())
  1602. {
  1603. psiEP_rec = qRecWestZDC[iHarm].Phi() / (iHarm + 1.);
  1604. if (psiEP_rec < 0.)
  1605. {
  1606. psiEP_rec += 2 * TMath::Pi() / (iHarm + 1.);
  1607. }
  1608. }
  1609. dPsi = 0.;
  1610. for (int iK = 0; iK < NShiftOrderMax * 3; iK++)
  1611. {
  1612. shiftCos = p_ZDC_Cos_West_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1613. shiftSin = p_ZDC_Sin_West_EP[iHarm][VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1614. dPsi += (1. / (iHarm + 1.)) * (2. / (iK + 1.)) *
  1615. (-1. * shiftSin * cos((iHarm + 1.) * (iK + 1.) * psiEP_rec) +
  1616. 1. * shiftCos * sin((iHarm + 1.) * (iK + 1.) * psiEP_rec));
  1617. }
  1618. psiEP_ZDC_West[iHarm] = psiEP_rec + dPsi;
  1619. if (psiEP_ZDC_West[iHarm] < 0.0)
  1620. psiEP_ZDC_West[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
  1621. if (psiEP_ZDC_West[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
  1622. psiEP_ZDC_West[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
  1623. }
  1624. // Track loop
  1625. for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
  1626. {
  1627. // Retrieve i-th femto track
  1628. StFemtoTrack *femtoTrack = dst->track(iTrk);
  1629. if (!femtoTrack)
  1630. continue;
  1631. // Must be a primary track
  1632. if (!femtoTrack->isPrimary())
  1633. continue;
  1634. //Track selection
  1635. if (!isGoodTrack(femtoTrack, energy, pVtx)) continue;
  1636. fUsedTracks++;
  1637. //Fill recentered Q-vectors for all TPC (FULL)
  1638. weight = GetWeight(femtoTrack);
  1639. QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1640. QvTMP.SetY(p_Qx_FULL_EP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1641. fQv2FullEP += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
  1642. QvTMP.Set(0., 0.);
  1643. QvTMP.SetX(p_Qx_FULL_EP[1][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1644. QvTMP.SetY(p_Qy_FULL_EP[1][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1645. fQv3FullEP += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
  1646. QvTMP.Set(0., 0.);
  1647. QvTMP.SetX(p_Qx_FULL_SP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1648. QvTMP.SetY(p_Qy_FULL_SP[0][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1649. fQv2FullSP += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
  1650. QvTMP.Set(0., 0.);
  1651. QvTMP.SetX(p_Qx_FULL_SP[1][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1652. QvTMP.SetY(p_Qy_FULL_SP[1][VtxSign][(Int_t) event->cent16()][(Double_t) event->runId()]);
  1653. fQv3FullSP += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
  1654. QvTMP.Set(0., 0.);
  1655. fQcountFull++;
  1656. if (femtoTrack->eta() >= 0.)
  1657. {
  1658. fQcountFullWest++;
  1659. } else
  1660. {
  1661. fQcountFullEast++;
  1662. }
  1663. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  1664. {
  1665. //fill Qvectors from EAST arm
  1666. if (femtoTrack->eta() > -1. * cutEta && femtoTrack->eta() < 0.)//cutEtaGap[iGap])
  1667. {
  1668. QvTMP.SetX(p_Qx_EAST_EP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1669. QvTMP.SetY(p_Qy_EAST_EP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1670. fQv2EastEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
  1671. QvTMP.Set(0., 0.);
  1672. QvTMP.SetX(p_Qx_EAST_EP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1673. QvTMP.SetY(p_Qy_EAST_EP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1674. fQv3EastEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
  1675. QvTMP.Set(0., 0.);
  1676. QvTMP.SetX(p_Qx_EAST_SP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1677. QvTMP.SetY(p_Qy_EAST_SP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1678. fQv2EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
  1679. QvTMP.Set(0., 0.);
  1680. QvTMP.SetX(p_Qx_EAST_SP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1681. QvTMP.SetY(p_Qy_EAST_SP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1682. fQv3EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
  1683. QvTMP.Set(0., 0.);
  1684. fQcountEast[iGap]++;
  1685. }
  1686. //fill Qvectors from WEST arm
  1687. if (femtoTrack->eta() < 1. * cutEta && femtoTrack->eta() > 0.)//cutEtaGap[iGap])
  1688. {
  1689. QvTMP.SetX(p_Qx_WEST_EP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1690. QvTMP.SetY(p_Qy_WEST_EP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1691. fQv2WestEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
  1692. QvTMP.Set(0., 0.);
  1693. QvTMP.SetX(p_Qx_WEST_EP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1694. QvTMP.SetY(p_Qy_WEST_EP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1695. fQv3WestEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
  1696. QvTMP.Set(0., 0.);
  1697. QvTMP.SetX(p_Qx_WEST_SP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1698. QvTMP.SetY(p_Qy_WEST_SP[0][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1699. fQv2WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
  1700. QvTMP.Set(0., 0.);
  1701. QvTMP.SetX(p_Qx_WEST_SP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1702. QvTMP.SetY(p_Qy_WEST_SP[1][VtxSign][(Int_t) event->cent16()][iGap][(Double_t) event->runId()]);
  1703. fQv3WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
  1704. QvTMP.Set(0., 0.);
  1705. fQcountWest[iGap]++;
  1706. }
  1707. }
  1708. // Check if track has TOF signal
  1709. if (femtoTrack->isTofTrack())
  1710. {
  1711. } //if( isTofTrack() )
  1712. } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
  1713. //Getting PsiEP with recentering + shift correction
  1714. Psi_ReCenter = 0.;
  1715. Double_t mCos, mSin;
  1716. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  1717. {
  1718. Psi2FullEP = 0.;
  1719. Psi3FullEP = 0.;
  1720. //FULL 2-nd order
  1721. Psi_ReCenter = TMath::ATan2(fQv2FullEP.Y(), fQv2FullEP.X()) / 2.;
  1722. dPsi = 0.;
  1723. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1724. {
  1725. mCos = p_Cos2_FULL_EP[VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1726. mSin = p_Cos2_FULL_EP[VtxSign][iK][(Int_t) event->cent16()][(Double_t) event->runId()];
  1727. dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) +
  1728. mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
  1729. }
  1730. Psi2FullEP = AngleShift(Psi_ReCenter + dPsi, 2.);
  1731. }
  1732. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  1733. {
  1734. if (fQcountEast[iGap] > cutNcountQvGap &&
  1735. fQcountWest[iGap] > cutNcountQvGap)
  1736. {
  1737. Psi2EastEP[iGap] = 0.;
  1738. Psi3EastEP[iGap] = 0.;
  1739. Psi2WestEP[iGap] = 0.;
  1740. Psi3WestEP[iGap] = 0.;
  1741. //EAST 2-nd order
  1742. Psi_ReCenter = TMath::ATan2(fQv2EastEP[iGap].Y(), fQv2EastEP[iGap].X()) / 2.;
  1743. dPsi = 0.;
  1744. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1745. {
  1746. mCos = p_Cos2_EAST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1747. mSin = p_Cos2_EAST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1748. dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin *
  1749. TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) +
  1750. mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
  1751. }
  1752. Psi2EastEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
  1753. //WEST 2-nd order
  1754. Psi_ReCenter = TMath::ATan2(fQv2WestEP[iGap].Y(), fQv2WestEP[iGap].X()) / 2.;
  1755. dPsi = 0.;
  1756. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1757. {
  1758. mCos = p_Cos2_WEST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1759. mSin = p_Cos2_WEST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1760. dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin *
  1761. TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) +
  1762. mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
  1763. }
  1764. Psi2WestEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
  1765. //EAST 3-nd order
  1766. Psi_ReCenter = TMath::ATan2(fQv3EastEP[iGap].Y(), fQv3EastEP[iGap].X()) / 3.;
  1767. dPsi = 0.;
  1768. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1769. {
  1770. mCos = p_Cos3_EAST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1771. mSin = p_Cos3_EAST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1772. dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin *
  1773. TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) +
  1774. mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
  1775. }
  1776. Psi3EastEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
  1777. //WEST 3-nd order
  1778. Psi_ReCenter = TMath::ATan2(fQv3WestEP[iGap].Y(), fQv3WestEP[iGap].X()) / 3.;
  1779. dPsi = 0.;
  1780. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1781. {
  1782. mCos = p_Cos3_WEST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1783. mSin = p_Cos3_WEST_EP[VtxSign][iK][(Int_t) event->cent16()][iGap][(Double_t) event->runId()];
  1784. dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin *
  1785. TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) +
  1786. mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
  1787. }
  1788. Psi3WestEP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
  1789. }
  1790. }
  1791. //FLOW CALCULATIONS
  1792. // Second Track loop
  1793. for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
  1794. {
  1795. // Retrieve i-th femto track
  1796. StFemtoTrack *femtoTrack = dst->track(iTrk);
  1797. if (!femtoTrack)
  1798. continue;
  1799. // Must be a primary track
  1800. if (!femtoTrack->isPrimary())
  1801. continue;
  1802. //Track selection
  1803. if (isGoodTrackFlowPID(femtoTrack, energy, pVtx))
  1804. //&& isGoodPID(femtoTrack))
  1805. {
  1806. weight = GetWeight(femtoTrack);
  1807. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  1808. {
  1809. res2EPsqr = p_res2_EP[iGap]->GetBinContent(p_res2_EP[iGap]->FindBin(event->cent16()));
  1810. res2EP = TMath::Sqrt(res2EPsqr);
  1811. res3EPsqr = p_res3_EP[iGap]->GetBinContent(p_res3_EP[iGap]->FindBin(event->cent16()));
  1812. res3EP = TMath::Sqrt(res3EPsqr);
  1813. //Loop over particle species (pi,K,p)
  1814. i_part = -1;
  1815. //TPC-only
  1816. if (!femtoTrack->isTofTrack())
  1817. {
  1818. //pion id
  1819. if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.6 &&
  1820. TMath::Abs(femtoTrack->nSigmaPion()) < 2)
  1821. {
  1822. i_part = 0;
  1823. }
  1824. // kaon id
  1825. if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.5 &&
  1826. TMath::Abs(femtoTrack->nSigmaKaon()) < 2)
  1827. {
  1828. i_part = 1;
  1829. }
  1830. // proton id
  1831. if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 0.9 &&
  1832. TMath::Abs(femtoTrack->nSigmaProton()) < 2)
  1833. {
  1834. i_part = 2;
  1835. }
  1836. //pion id
  1837. if (femtoTrack->ptot() >= 0.6 && femtoTrack->ptot() < 0.7 &&
  1838. TMath::Abs(femtoTrack->nSigmaPion()) < 2 &&
  1839. TMath::Abs(femtoTrack->nSigmaKaon()) > 2)
  1840. {
  1841. i_part = 0;
  1842. }
  1843. // kaon id
  1844. if (femtoTrack->ptot() >= 0.5 && femtoTrack->ptot() < 0.7 &&
  1845. TMath::Abs(femtoTrack->nSigmaKaon()) < 2 &&
  1846. TMath::Abs(femtoTrack->nSigmaPion()) > 3)
  1847. {
  1848. i_part = 1;
  1849. }
  1850. // proton id
  1851. if (femtoTrack->ptot() >= 0.9 && femtoTrack->ptot() < 1.2 &&
  1852. TMath::Abs(femtoTrack->nSigmaProton()) < 2 &&
  1853. TMath::Abs(femtoTrack->nSigmaPion()) > 3)
  1854. {
  1855. i_part = 2;
  1856. }
  1857. }
  1858. //TPC+TOF
  1859. if (isGoodPID(femtoTrack))
  1860. {
  1861. // pion id - asymmetry cut
  1862. if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
  1863. TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
  1864. femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)// &&
  1865. {
  1866. i_part = 0;
  1867. }
  1868. // kaon id - asymmetry cut
  1869. if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
  1870. TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
  1871. femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)// &&
  1872. {
  1873. i_part = 1;
  1874. }
  1875. // proton id
  1876. if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
  1877. TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
  1878. femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)// &&
  1879. {
  1880. i_part = 2;
  1881. }
  1882. }
  1883. //}
  1884. if (i_part == -1) continue;
  1885. //EAST
  1886. //if (track->eta() > -1.*cutEta && track->eta() < 0.)
  1887. if (femtoTrack->eta() > -1. * cutEta && femtoTrack->eta() < -1. * cutEtaGap.at(iGap))
  1888. {
  1889. flow2EP = TMath::Cos(2. * (femtoTrack->phi() - Psi2WestEP[iGap])) / res2EP;
  1890. flow3EP = TMath::Cos(3. * (femtoTrack->phi() - Psi3WestEP[iGap])) / res3EP;
  1891. p_v2run_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
  1892. p_v3run_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP, weight);
  1893. }
  1894. //WEST
  1895. if (femtoTrack->eta() < 1. * cutEta && femtoTrack->eta() > 1. * cutEtaGap.at(iGap))
  1896. {
  1897. flow2EP = TMath::Cos(2. * (femtoTrack->phi() - Psi2EastEP[iGap])) / res2EP;
  1898. flow3EP = TMath::Cos(3. * (femtoTrack->phi() - Psi3EastEP[iGap])) / res3EP;
  1899. p_v2run_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
  1900. p_v3run_EP[i_part][iGap][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP, weight);
  1901. }
  1902. if (iGap == 0)
  1903. {
  1904. // Flow with BBC EP
  1905. res2EPBBC = hResBBCFULL[1]->GetBinContent(hResBBCFULL[1]->FindBin(event->cent16()));
  1906. flow2EPBBC = TMath::Cos(2. * (femtoTrack->phi() - psiEP_BBC_Full[1])) / res2EPBBC;
  1907. p_v2run_BBC[i_part][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EPBBC, weight);
  1908. // Flow with ZDC EP
  1909. // Full
  1910. res2EPZDC = hResZDCFULL[1]->GetBinContent(hResZDCFULL[1]->FindBin(event->cent16()));
  1911. flow2EPZDC = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_Full[0])) / res2EPZDC;
  1912. p_v2run_ZDC[i_part][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EPZDC, weight);
  1913. // East
  1914. res2EPZDC = hResZDCHalf[1]->GetBinContent(hResZDCHalf[1]->FindBin(event->cent16()));
  1915. flow2EPZDC = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_East[0])) / res2EPZDC;
  1916. p_v2run_ZDC_East[i_part][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EPZDC,
  1917. weight);
  1918. // West
  1919. res2EPZDC = hResZDCHalf[1]->GetBinContent(hResZDCHalf[1]->FindBin(event->cent16()));
  1920. flow2EPZDC = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_West[0])) / res2EPZDC;
  1921. p_v2run_ZDC_West[i_part][(Int_t) event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EPZDC,
  1922. weight);
  1923. }
  1924. }
  1925. }
  1926. }
  1927. } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
  1928. femtoReader->Finish();
  1929. TFile *output = new TFile(outFile, "recreate");
  1930. for (int iPID = 0; iPID < 3; iPID++)
  1931. {
  1932. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  1933. {
  1934. for (int iCent = 0; iCent < 16; iCent++)
  1935. {
  1936. p_v2run_EP[iPID][iGap][iCent]->Write();
  1937. p_v3run_EP[iPID][iGap][iCent]->Write();
  1938. }
  1939. }
  1940. }
  1941. for (int iPID = 0; iPID < 3; iPID++)
  1942. {
  1943. for (int iCent = 0; iCent < 16; iCent++)
  1944. {
  1945. p_v2run_BBC[iPID][iCent]->Write();
  1946. }
  1947. }
  1948. for (int iPID = 0; iPID < 3; iPID++)
  1949. {
  1950. for (int iCent = 0; iCent < 16; iCent++)
  1951. {
  1952. p_v2run_ZDC[iPID][iCent]->Write();
  1953. p_v2run_ZDC_East[iPID][iCent]->Write();
  1954. p_v2run_ZDC_West[iPID][iCent]->Write();
  1955. }
  1956. }
  1957. output->Close();
  1958. std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
  1959. << std::endl;
  1960. std::cout << "Good events: " << goodEventCounter << std::endl;
  1961. timer.Stop();
  1962. timer.Print();
  1963. }
  1964. Bool_t isGoodEvent(StFemtoEvent *const &event)
  1965. {
  1966. if (!event) return false;
  1967. if (event == nullptr) return false;
  1968. if (event->primaryVertex().Perp() > cutVtxR) return false;
  1969. if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy)) return false;
  1970. if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz) return false;
  1971. return true;
  1972. }
  1973. Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
  1974. {
  1975. if (!track) return false;
  1976. // if (!track->isPrimary()) return false;
  1977. if (track->nHitsFit() < cutNhits) return false;
  1978. if (track->nHitsPoss() <= cutNhitsPoss) return false;
  1979. if ((Double_t) track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio) return false;
  1980. if (TMath::Abs(track->eta()) >= cutEta) return false;
  1981. if (track->pt() <= cutPtMin.at(_energy)) return false;
  1982. if (track->pt() > cutPtMax) return false;
  1983. if (track->ptot() > cutPMax) return false;
  1984. if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy)) return false;
  1985. return true;
  1986. }
  1987. Bool_t isGoodTrackFlowPID(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
  1988. {
  1989. if (!track) return false;
  1990. // if (!track->isPrimary()) return false;
  1991. if (track->nHitsFit() < cutNhits) return false;
  1992. if (track->nHitsPoss() <= cutNhitsPoss) return false;
  1993. if ((Double_t) track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio) return false;
  1994. if (TMath::Abs(track->eta()) >= cutEta) return false;
  1995. if (track->pt() <= cutPtMin.at(_energy)) return false;
  1996. //if (track->pt() > 4.15) return false;
  1997. if (track->ptot() > cutPMax) return false;
  1998. if (track->gDCA(pVtx).Mag() >= cutDCA_PID) return false;
  1999. return true;
  2000. }
  2001. Bool_t isGoodTrackFlowPIDsys_Eta(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount)
  2002. {
  2003. if (!track) return false;
  2004. // if (!track->isPrimary()) return false;
  2005. if (track->nHitsFit() < cutNhits) return false;
  2006. if (track->nHitsPoss() <= cutNhitsPoss) return false;
  2007. if ((Double_t) track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio) return false;
  2008. if (TMath::Abs(track->eta()) >= cutEtasys.at(sysCount)) return false;
  2009. if (track->pt() <= cutPtMin.at(_energy)) return false;
  2010. //if (track->pt() > 4.15) return false;
  2011. if (track->ptot() > cutPMax) return false;
  2012. if (track->gDCA(pVtx).Mag() >= cutDCA_PID) return false;
  2013. return true;
  2014. }
  2015. Bool_t isGoodTrackFlowPIDsys_Nhits(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount)
  2016. {
  2017. if (!track) return false;
  2018. // if (!track->isPrimary()) return false;
  2019. if (track->nHitsFit() < cutNhitssys.at(sysCount)) return false;
  2020. if (track->nHitsPoss() <= cutNhitsPoss) return false;
  2021. if (TMath::Abs(track->eta()) >= cutEta) return false;
  2022. if (track->pt() <= cutPtMin.at(_energy)) return false;
  2023. //if (track->pt() > 4.15) return false;
  2024. if (track->ptot() > cutPMax) return false;
  2025. if (track->gDCA(pVtx).Mag() >= cutDCA_PID) return false;
  2026. return true;
  2027. }
  2028. Bool_t isGoodTrackFlowPIDsys_DCA(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx, Int_t sysCount)
  2029. {
  2030. if (!track) return false;
  2031. // if (!track->isPrimary()) return false;
  2032. if (track->nHitsFit() < cutNhits) return false;
  2033. if (track->nHitsPoss() <= cutNhitsPoss) return false;
  2034. if ((Double_t) track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio) return false;
  2035. if (TMath::Abs(track->eta()) >= cutEta) return false;
  2036. if (track->pt() <= cutPtMin.at(_energy)) return false;
  2037. //if (track->pt() > 4.15) return false;
  2038. if (track->ptot() > cutPMax) return false;
  2039. if (track->gDCA(pVtx).Mag() >= cutDCA_PIDsys.at(sysCount)) return false;
  2040. return true;
  2041. }
  2042. Bool_t isGoodPID(StFemtoTrack *const &track)
  2043. {
  2044. if (!track->isTofTrack()) return false;
  2045. if (track->massSqr() < cutMass2Min) return false;
  2046. //NhitsDedx cut applied in StFemtoDst?
  2047. // ToFYLocal cut applied in StFemtoDst?
  2048. return true;
  2049. }
  2050. Double_t GetWeight(StFemtoTrack *const &track)
  2051. {
  2052. Double_t weight;
  2053. if (track->pt() <= cutPtWeightEP)
  2054. {
  2055. weight = track->pt();
  2056. } else
  2057. {
  2058. weight = cutPtWeightEP;
  2059. }
  2060. return weight;
  2061. }
  2062. TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
  2063. {
  2064. TVector2 qv(0., 0.);
  2065. qv.Set(TMath::Cos(_harm * track->phi()), TMath::Sin(_harm * track->phi()));
  2066. return qv;
  2067. }
  2068. Double_t AngleShift(Double_t Psi, Double_t order)
  2069. {
  2070. Double_t PsiCorr = Psi;
  2071. if (PsiCorr > TMath::Pi() / order)
  2072. {
  2073. PsiCorr = Psi - 2. * TMath::Pi() / order;
  2074. }
  2075. if (PsiCorr < -1. * TMath::Pi() / order)
  2076. {
  2077. PsiCorr = Psi + 2. * TMath::Pi() / order;
  2078. }
  2079. return PsiCorr;
  2080. }
  2081. Int_t GetVzBin(Double_t vtxZ)
  2082. {
  2083. for (Int_t i = 0; i < fArraySize; i++)
  2084. {
  2085. if (vtxZ >= VtxBins[i] && vtxZ < VtxBins[i + 1]) return i;
  2086. }
  2087. return -1;
  2088. }
  2089. //Int_t DoPID(StFemtoTrack *const &track)
  2090. //{
  2091. // if (track->isTofTrack())
  2092. // {
  2093. // if (track->ptot() >= 0.2 && track->ptot() < 3.4 &&
  2094. // TMath::Abs(track->nSigmaPion()) < 3. &&
  2095. // track->massSqr() > pionMassSqr - NpidFactor*sigMsqrPion->Eval(femtoTrack->pt()) &&
  2096. // femtoTrack->massSqr() < pionMassSqr + NpidFactor*sigMsqrPion->Eval(femtoTrack->pt()))
  2097. // }
  2098. //}
  2099. Double_t GetSigmWidth(StFemtoTrack *const &track, Int_t PID)
  2100. {
  2101. Double_t pol0, pol1, pol2, pol3;
  2102. Double_t pt = track->pt();
  2103. if (PID != 0 && PID != 1 && PID != 2) return 0.;
  2104. if (PID == 0) // pion
  2105. {
  2106. pol0 = -0.036277878;
  2107. pol1 = 0.10230040;
  2108. pol2 = -0.069674611;
  2109. pol3 = 0.020434016;
  2110. }
  2111. if (PID == 1) // kaon
  2112. {
  2113. pol0 = 0.011004116;
  2114. pol1 = -0.022878517;
  2115. pol2 = 0.038562140;
  2116. pol3 = -0.0060050387;
  2117. }
  2118. if (PID == 2) // proton
  2119. {
  2120. pol0 = 0.035215709;
  2121. pol1 = -0.038627269;
  2122. pol2 = 0.046936935;
  2123. pol3 = -0.0063598115;
  2124. }
  2125. return (pol0 + pol1 * pt + pol2 * pt * pt + pol3 * pt * pt * pt);
  2126. }
  2127. Double_t GetMeanNsig(StFemtoTrack *const &track, const Int_t PID)
  2128. {
  2129. Double_t pt, par0, par1, par2, par3;
  2130. if (PID == 0) return 0.;
  2131. if (PID == 1)
  2132. {
  2133. par0 = 5.84263e+00;
  2134. par1 = -6.65912e+00;
  2135. par2 = 8.32204e-01;
  2136. par3 = 8.28440e-01;
  2137. }
  2138. if (PID == 2)
  2139. {
  2140. par0 = 1.09685e+01;
  2141. par1 = -8.61499e+00;
  2142. par2 = 8.92938e-01;
  2143. par3 = 8.08397e-01;
  2144. }
  2145. pt = track->pt();
  2146. if (pt != 0)
  2147. {
  2148. return (par0 * TMath::Power(1. / pt, par2) + par1 + par3 * pt);
  2149. } else
  2150. {
  2151. return 0.;
  2152. }
  2153. }
  2154. Double_t GetWidthM2(StFemtoTrack *const &track, const Int_t PID)
  2155. {
  2156. Double_t pt, pol0, pol1, pol2;
  2157. if (PID != 0) return 0.;
  2158. pol0 = 8.74975e-04;
  2159. pol1 = -1.62659e-03;
  2160. pol2 = 2.89828e-02;
  2161. pt = track->pt();
  2162. return (pol0 + pol1 * pt + pol2 * pt * pt);
  2163. }
  2164. Double_t GetScaleFactor(StFemtoTrack *const &track)
  2165. {
  2166. Double_t wM2pion = GetWidthM2(track, 0);
  2167. Double_t exScale = 1. / 0.8273; // = 1. for BES energies
  2168. if (wM2pion == 0) return 1.;
  2169. else return (exScale / wM2pion);
  2170. }
  2171. Double_t GetRotAngle(StFemtoTrack *const &track)
  2172. {
  2173. Double_t scaleFactor = GetScaleFactor(track);
  2174. Double_t MeanKaonNsig = GetMeanNsig(track, 1);
  2175. Double_t new_x = (MeanKaonNsig - 0.) / scaleFactor;
  2176. Double_t new_y = (kaonMassSqr - pionMassSqr);
  2177. return (-TMath::ATan2(new_y, new_x));
  2178. }
  2179. std::pair<Double_t, Double_t> GetNewXY(StFemtoTrack *const &track)
  2180. {
  2181. std::pair<Double_t, Double_t> coord;
  2182. Double_t scaleFactor = GetScaleFactor(track);
  2183. Double_t new_x = (track->nSigmaPion() - 0.) / scaleFactor;
  2184. Double_t new_y = (track->massSqr() - pionMassSqr);
  2185. Double_t rotAngle = GetRotAngle(track);
  2186. coord.first = new_x * TMath::Cos(rotAngle) - new_y * TMath::Sin(rotAngle);
  2187. coord.second = new_x * TMath::Sin(rotAngle) + new_y * TMath::Cos(rotAngle);
  2188. return coord;
  2189. }
  2190. Bool_t isMesonXY(Double_t x, Double_t y)
  2191. // Used to exclude proton peak from new XY coord distribution
  2192. {
  2193. Double_t pol0, pol1;
  2194. pol0 = -0.92857143;
  2195. pol1 = 1.4285714;
  2196. // Not a proton
  2197. if (y > (pol0 + pol1 * x)) return true;
  2198. // Proton
  2199. if (y <= (pol0 + pol1 * x)) return false;
  2200. // just in case...
  2201. return false;
  2202. }
  2203. Double_t GetSigmNewX(StFemtoTrack *const &track, Int_t PID)
  2204. {
  2205. if (PID != 0 && PID != 1 && PID != 2) return 0.;
  2206. Double_t pol0, pol1, pol2, pol3;
  2207. Double_t pt = track->pt();
  2208. if (PID == 0) // pion
  2209. {
  2210. pol0 = 0.00061920122;
  2211. pol1 = 0.00064431355;
  2212. pol2 = 0.013672155;
  2213. pol3 = 0.0012772833;
  2214. }
  2215. if (PID == 1) // kaon
  2216. {
  2217. pol0 = 0.0057203659;
  2218. pol1 = 0.0013767144;
  2219. pol2 = 0.0063761902;
  2220. pol3 = 0.0073111853;
  2221. }
  2222. return (pol0 + pol1 * pt + pol2 * pt * pt + pol3 * pt * pt * pt);
  2223. }
  2224. Double_t GetSigmNewY(StFemtoTrack *const &track, Int_t PID)
  2225. {
  2226. if (PID != 0 && PID != 1 && PID != 2) return 0.;
  2227. Double_t pol0, pol1, pol2, pol3;
  2228. Double_t pt = track->pt();
  2229. if (PID == 0) // pion
  2230. {
  2231. pol0 = 0.0017345792;
  2232. pol1 = -0.0059824003;
  2233. pol2 = 0.029537517;
  2234. pol3 = -0.0037990834;
  2235. }
  2236. if (PID == 1) // kaon
  2237. {
  2238. pol0 = 0.0041058086;
  2239. pol1 = -0.0041708451;
  2240. pol2 = 0.026109827;
  2241. pol3 = -0.0018772891;
  2242. }
  2243. return (pol0 + pol1 * pt + pol2 * pt * pt + pol3 * pt * pt * pt);
  2244. }
  2245. //--------------------------------------------------------------------------------------------------------------------------//
  2246. //this function simply connects the gain values read in to the BBC azimuthal distribution
  2247. //since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
  2248. //I prefer assigning the angle between the tiles thus "greying out" the adcs.
  2249. //Others have assigned all of the adc to one (exclusive) or the the other.
  2250. Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId)
  2251. {
  2252. //float GetPhiInBBC(int eastWest, int bbcN) { //tileId=0 to 23
  2253. const float Pi = TMath::Pi();
  2254. const float Phi_div = Pi / 6;
  2255. float bbc_phi = Phi_div;
  2256. switch (tileId)
  2257. {
  2258. case 0:
  2259. bbc_phi = 3. * Phi_div;
  2260. break;
  2261. case 1:
  2262. bbc_phi = Phi_div;
  2263. break;
  2264. case 2:
  2265. bbc_phi = -1. * Phi_div;
  2266. break;
  2267. case 3:
  2268. bbc_phi = -3. * Phi_div;
  2269. break;
  2270. case 4:
  2271. bbc_phi = -5. * Phi_div;
  2272. break;
  2273. case 5:
  2274. bbc_phi = 5. * Phi_div;
  2275. break;
  2276. //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
  2277. case 6:
  2278. bbc_phi = 3. * Phi_div;
  2279. break;
  2280. case 7:
  2281. bbc_phi = 3. * Phi_div;
  2282. break;
  2283. case 8:
  2284. bbc_phi = Phi_div;
  2285. break;
  2286. case 9:
  2287. bbc_phi = 0.;
  2288. break;
  2289. case 10:
  2290. bbc_phi = -1. * Phi_div;
  2291. break;
  2292. //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div; //tiles 13 and 15 are gained together
  2293. case 11:
  2294. bbc_phi = -3. * Phi_div;
  2295. break;
  2296. case 12:
  2297. bbc_phi = -3. * Phi_div;
  2298. break;
  2299. case 13:
  2300. bbc_phi = -5. * Phi_div;
  2301. break;
  2302. case 14:
  2303. bbc_phi = Pi;
  2304. break;
  2305. case 15:
  2306. bbc_phi = 5. * Phi_div;
  2307. break;
  2308. }
  2309. //if we're looking at the east BBC we need to flip around x in the STAR coordinates,
  2310. //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the
  2311. if (0 == eastWest)
  2312. {
  2313. if (bbc_phi > -0.001)
  2314. { //this is not a >= since we are talking about finite adcs -- not to important
  2315. bbc_phi = Pi - bbc_phi;
  2316. } else
  2317. {
  2318. bbc_phi = -Pi - bbc_phi;
  2319. }
  2320. }
  2321. if (bbc_phi < 0.0)
  2322. bbc_phi += 2. * Pi;
  2323. if (bbc_phi > 2. * Pi)
  2324. bbc_phi -= 2. * Pi;
  2325. return bbc_phi;
  2326. }
  2327. Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip)
  2328. // Get position of each slat;strip starts from 0
  2329. {
  2330. Double_t zdcsmd_x[7] = {0.5, 2, 3.5, 5, 6.5, 8, 9.5};
  2331. Double_t zdcsmd_y[8] = {1.25, 3.25, 5.25, 7.25, 9.25, 11.25, 13.25, 15.25};
  2332. Double_t mZDCSMDCenterex = 4.72466;
  2333. Double_t mZDCSMDCenterey = 5.53629;
  2334. Double_t mZDCSMDCenterwx = 4.39604;
  2335. Double_t mZDCSMDCenterwy = 5.19968;
  2336. if (eastwest == 0 && verthori == 0)
  2337. return zdcsmd_x[strip] - mZDCSMDCenterex;
  2338. if (eastwest == 1 && verthori == 0)
  2339. return mZDCSMDCenterwx - zdcsmd_x[strip];
  2340. if (eastwest == 0 && verthori == 1)
  2341. return zdcsmd_y[strip] / sqrt(2.) - mZDCSMDCenterey;
  2342. if (eastwest == 1 && verthori == 1)
  2343. return zdcsmd_y[strip] / sqrt(2.) - mZDCSMDCenterwy;
  2344. return -999.;
  2345. }