FemtoDstAnalyzer_QA.C 119 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483248424852486248724882489249024912492249324942495249624972498249925002501250225032504250525062507250825092510251125122513251425152516251725182519252025212522252325242525252625272528252925302531253225332534253525362537253825392540254125422543254425452546254725482549255025512552255325542555255625572558255925602561256225632564256525662567256825692570257125722573257425752576257725782579258025812582258325842585258625872588258925902591259225932594259525962597259825992600260126022603260426052606260726082609261026112612261326142615261626172618261926202621262226232624262526262627262826292630263126322633263426352636263726382639264026412642264326442645264626472648264926502651265226532654265526562657265826592660266126622663266426652666266726682669267026712672267326742675267626772678267926802681268226832684268526862687268826892690269126922693269426952696269726982699270027012702270327042705270627072708270927102711271227132714271527162717
  1. /**
  2. * \brief Example of how to read a file (list of files) using StFemtoEvent classes
  3. *
  4. * RunFemtoDstAnalyzer.C is an example of reading FemtoDst format.
  5. * One can use either FemtoDst file or a list of femtoDst files (inFile.lis or
  6. * inFile.list) as an input, and preform physics analysis
  7. *
  8. * \author Grigory Nigmatkulov
  9. * \date May 29, 2018
  10. */
  11. // This is needed for calling standalone classes
  12. #define _VANILLA_ROOT_
  13. // C++ headers
  14. #include <iostream>
  15. #include <vector>
  16. #include <map>
  17. // ROOT headers
  18. #include "TROOT.h"
  19. #include "TFile.h"
  20. #include "TChain.h"
  21. #include "TVector2.h"
  22. #include "TVector3.h"
  23. #include "TTree.h"
  24. #include "TSystem.h"
  25. #include "TH1.h"
  26. #include "TH2.h"
  27. #include "TMath.h"
  28. #include "TProfile2D.h"
  29. #include "TStopwatch.h"
  30. #include "TRandom3.h"
  31. #include "TVector2.h"
  32. // FemtoDst headers
  33. #include "StFemtoDstReader.h"
  34. #include "StFemtoDst.h"
  35. #include "StFemtoEvent.h"
  36. #include "StFemtoTrack.h"
  37. // Load libraries (for ROOT_VERSTION_CODE >= 393215)
  38. #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
  39. R__LOAD_LIBRARY(/ mnt / pool / rhic / 4 / parfenovpeter / STAR / build / libStFemtoDst.so)
  40. #endif
  41. #include "Constants.h"
  42. // inFile - is a name of name.FemtoDst.root file or a name
  43. // of a name.lis(t) files, that contains a list of
  44. // name1.FemtoDst.root, name2.FemtoDst.root, ... files
  45. //Used functions (see them below)
  46. Bool_t isGoodEvent(StFemtoEvent *const &event);
  47. Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
  48. Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
  49. Double_t GetWeight(StFemtoTrack *const &track);
  50. TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
  51. Double_t AngleShift(Double_t Psi, Double_t order);
  52. Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);
  53. Int_t GetVzBin(Double_t vtxZ);
  54. Double_t GetEPAngle(TVector2 qvect, Double_t order);
  55. Double_t GetEPAngle(Double_t angle, Double_t order);
  56. Bool_t isGoodPID(StFemtoTrack *const &track);
  57. //_________________
  58. void FemtoDstAnalyzer_QA(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
  59. const Char_t *outFile = "./oShiftTest.root",
  60. const Char_t *reCentFile = "",
  61. const Char_t *shiftFile = "",
  62. const Char_t *gainBBC = "",
  63. const Char_t *reCentFileBBC = "",
  64. const Char_t *shiftFileBBC = "")
  65. {
  66. std::cout << "Hi! Lets do some physics, Master!" << std::endl;
  67. TStopwatch timer;
  68. timer.Start();
  69. #if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
  70. gSystem->Load("../libStFemtoDst.so");
  71. #endif
  72. StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
  73. femtoReader->Init();
  74. // This is a way if you want to spead up IO
  75. std::cout << "Explicit read status for some branches" << std::endl;
  76. femtoReader->SetStatus("*", 0);
  77. femtoReader->SetStatus("Event", 1);
  78. femtoReader->SetStatus("Track", 1);
  79. std::cout << "Status has been set" << std::endl;
  80. std::cout << "Now I know what to read, Master!" << std::endl;
  81. if (!femtoReader->chain())
  82. {
  83. std::cout << "No chain has been found." << std::endl;
  84. }
  85. Long64_t eventsInTree = femtoReader->tree()->GetEntries();
  86. std::cout << "eventsInTree: " << eventsInTree << std::endl;
  87. Long64_t events2read = femtoReader->chain()->GetEntries();
  88. std::cout << "Number of events to read: " << events2read << std::endl;
  89. //Read ReCentering <Q> vectors
  90. TFile *fiReCent = new TFile(reCentFile, "read");
  91. std::map<Double_t, Double_t> p_Qx_FULL_EP[fNharmonics][fArraySize][fNCent];
  92. std::map<Double_t, Double_t> p_Qy_FULL_EP[fNharmonics][fArraySize][fNCent];
  93. std::map<Double_t, Double_t> p_Qx_FULL_SP[fNharmonics][fArraySize][fNCent];
  94. std::map<Double_t, Double_t> p_Qy_FULL_SP[fNharmonics][fArraySize][fNCent];
  95. std::map<Double_t, Double_t> p_Qx_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  96. std::map<Double_t, Double_t> p_Qy_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  97. std::map<Double_t, Double_t> p_Qx_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  98. std::map<Double_t, Double_t> p_Qy_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  99. std::map<Double_t, Double_t> p_Qx_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  100. std::map<Double_t, Double_t> p_Qy_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  101. std::map<Double_t, Double_t> p_Qx_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  102. std::map<Double_t, Double_t> p_Qy_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
  103. TProfile2D *prof;
  104. double binCont;
  105. std::cout << "Reading file for TPC recentering..." << std::endl;
  106. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  107. {
  108. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  109. {
  110. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
  111. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_EP%i", iHarm, iVtx));
  112. for (int i = 0; i < fNBins; i++)
  113. {
  114. for (int iCent = 0; iCent < fNCent; iCent++)
  115. {
  116. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  117. if (binCont == 0.)
  118. continue;
  119. p_Qx_FULL_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  120. }
  121. }
  122. delete prof;
  123. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_EP%i", iHarm, iVtx));
  124. for (int i = 0; i < fNBins; i++)
  125. {
  126. for (int iCent = 0; iCent < fNCent; iCent++)
  127. {
  128. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  129. if (binCont == 0.)
  130. continue;
  131. p_Qy_FULL_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  132. }
  133. }
  134. delete prof;
  135. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_SP%i", iHarm, iVtx));
  136. for (int i = 0; i < fNBins; i++)
  137. {
  138. for (int iCent = 0; iCent < fNCent; iCent++)
  139. {
  140. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  141. if (binCont == 0.)
  142. continue;
  143. p_Qx_FULL_SP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  144. }
  145. }
  146. delete prof;
  147. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_SP%i", iHarm, iVtx));
  148. for (int i = 0; i < fNBins; i++)
  149. {
  150. for (int iCent = 0; iCent < fNCent; iCent++)
  151. {
  152. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  153. if (binCont == 0.)
  154. continue;
  155. p_Qy_FULL_SP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  156. }
  157. }
  158. delete prof;
  159. for (int iEtaGap = 0; iEtaGap < NEtaGaps; iEtaGap++)
  160. {
  161. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  162. for (int i = 0; i < fNBins; i++)
  163. {
  164. for (int iCent = 0; iCent < fNCent; iCent++)
  165. {
  166. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  167. if (binCont == 0.)
  168. continue;
  169. p_Qx_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  170. }
  171. }
  172. delete prof;
  173. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  174. for (int i = 0; i < fNBins; i++)
  175. {
  176. for (int iCent = 0; iCent < fNCent; iCent++)
  177. {
  178. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  179. if (binCont == 0.)
  180. continue;
  181. p_Qy_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  182. }
  183. }
  184. delete prof;
  185. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  186. for (int i = 0; i < fNBins; i++)
  187. {
  188. for (int iCent = 0; iCent < fNCent; iCent++)
  189. {
  190. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  191. if (binCont == 0.)
  192. continue;
  193. p_Qx_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  194. }
  195. }
  196. delete prof;
  197. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  198. for (int i = 0; i < fNBins; i++)
  199. {
  200. for (int iCent = 0; iCent < fNCent; iCent++)
  201. {
  202. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  203. if (binCont == 0.)
  204. continue;
  205. p_Qy_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  206. }
  207. }
  208. delete prof;
  209. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  210. for (int i = 0; i < fNBins; i++)
  211. {
  212. for (int iCent = 0; iCent < fNCent; iCent++)
  213. {
  214. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  215. if (binCont == 0.)
  216. continue;
  217. p_Qx_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  218. }
  219. }
  220. delete prof;
  221. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
  222. for (int i = 0; i < fNBins; i++)
  223. {
  224. for (int iCent = 0; iCent < fNCent; iCent++)
  225. {
  226. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  227. if (binCont == 0.)
  228. continue;
  229. p_Qy_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  230. }
  231. }
  232. delete prof;
  233. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  234. for (int i = 0; i < fNBins; i++)
  235. {
  236. for (int iCent = 0; iCent < fNCent; iCent++)
  237. {
  238. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  239. if (binCont == 0.)
  240. continue;
  241. p_Qx_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  242. }
  243. }
  244. delete prof;
  245. prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
  246. for (int i = 0; i < fNBins; i++)
  247. {
  248. for (int iCent = 0; iCent < fNCent; iCent++)
  249. {
  250. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  251. if (binCont == 0.)
  252. continue;
  253. p_Qy_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
  254. }
  255. }
  256. delete prof;
  257. }
  258. }
  259. }
  260. //Read Shift coefficients
  261. TFile *fiShift = new TFile(shiftFile, "read");
  262. std::map<Double_t, Double_t> p_Cos2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  263. std::map<Double_t, Double_t> p_Sin2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  264. std::map<Double_t, Double_t> p_Cos2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  265. std::map<Double_t, Double_t> p_Sin2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  266. std::map<Double_t, Double_t> p_Cos2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  267. std::map<Double_t, Double_t> p_Sin2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  268. std::map<Double_t, Double_t> p_Cos2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  269. std::map<Double_t, Double_t> p_Sin2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  270. std::map<Double_t, Double_t> p_Cos2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  271. std::map<Double_t, Double_t> p_Sin2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  272. std::map<Double_t, Double_t> p_Cos2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  273. std::map<Double_t, Double_t> p_Sin2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  274. std::map<Double_t, Double_t> p_Cos3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  275. std::map<Double_t, Double_t> p_Sin3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
  276. std::map<Double_t, Double_t> p_Cos3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  277. std::map<Double_t, Double_t> p_Sin3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
  278. std::map<Double_t, Double_t> p_Cos3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  279. std::map<Double_t, Double_t> p_Sin3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  280. std::map<Double_t, Double_t> p_Cos3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  281. std::map<Double_t, Double_t> p_Sin3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  282. std::map<Double_t, Double_t> p_Cos3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  283. std::map<Double_t, Double_t> p_Sin3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  284. std::map<Double_t, Double_t> p_Cos3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  285. std::map<Double_t, Double_t> p_Sin3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
  286. std::cout << "Reading file for TPC flattening..." << std::endl;
  287. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  288. {
  289. for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
  290. {
  291. std::cout << "iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
  292. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
  293. for (int i = 0; i < fNBins; i++)
  294. {
  295. for (int iCent = 0; iCent < fNCent; iCent++)
  296. {
  297. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  298. if (binCont == 0.)
  299. continue;
  300. p_Cos2_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  301. }
  302. }
  303. delete prof;
  304. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
  305. for (int i = 0; i < fNBins; i++)
  306. {
  307. for (int iCent = 0; iCent < fNCent; iCent++)
  308. {
  309. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  310. if (binCont == 0.)
  311. continue;
  312. p_Sin2_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  313. }
  314. }
  315. delete prof;
  316. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
  317. for (int i = 0; i < fNBins; i++)
  318. {
  319. for (int iCent = 0; iCent < fNCent; iCent++)
  320. {
  321. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  322. if (binCont == 0.)
  323. continue;
  324. p_Cos2_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  325. }
  326. }
  327. delete prof;
  328. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
  329. for (int i = 0; i < fNBins; i++)
  330. {
  331. for (int iCent = 0; iCent < fNCent; iCent++)
  332. {
  333. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  334. if (binCont == 0.)
  335. continue;
  336. p_Sin2_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  337. }
  338. }
  339. delete prof;
  340. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
  341. for (int i = 0; i < fNBins; i++)
  342. {
  343. for (int iCent = 0; iCent < fNCent; iCent++)
  344. {
  345. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  346. if (binCont == 0.)
  347. continue;
  348. p_Cos3_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  349. }
  350. }
  351. delete prof;
  352. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
  353. for (int i = 0; i < fNBins; i++)
  354. {
  355. for (int iCent = 0; iCent < fNCent; iCent++)
  356. {
  357. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  358. if (binCont == 0.)
  359. continue;
  360. p_Sin3_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  361. }
  362. }
  363. delete prof;
  364. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
  365. for (int i = 0; i < fNBins; i++)
  366. {
  367. for (int iCent = 0; iCent < fNCent; iCent++)
  368. {
  369. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  370. if (binCont == 0.)
  371. continue;
  372. p_Cos3_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  373. }
  374. }
  375. delete prof;
  376. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
  377. for (int i = 0; i < fNBins; i++)
  378. {
  379. for (int iCent = 0; iCent < fNCent; iCent++)
  380. {
  381. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  382. if (binCont == 0.)
  383. continue;
  384. p_Sin3_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  385. }
  386. }
  387. delete prof;
  388. for (int iGap = 0; iGap < NEtaGaps; iGap++)
  389. {
  390. // for 2nd harmonics
  391. // EAST
  392. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  393. for (int i = 0; i < fNBins; i++)
  394. {
  395. for (int iCent = 0; iCent < fNCent; iCent++)
  396. {
  397. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  398. if (binCont == 0.)
  399. continue;
  400. p_Cos2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  401. }
  402. }
  403. delete prof;
  404. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  405. for (int i = 0; i < fNBins; i++)
  406. {
  407. for (int iCent = 0; iCent < fNCent; iCent++)
  408. {
  409. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  410. if (binCont == 0.)
  411. continue;
  412. p_Sin2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  413. }
  414. }
  415. delete prof;
  416. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  417. for (int i = 0; i < fNBins; i++)
  418. {
  419. for (int iCent = 0; iCent < fNCent; iCent++)
  420. {
  421. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  422. if (binCont == 0.)
  423. continue;
  424. p_Cos2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  425. }
  426. }
  427. delete prof;
  428. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  429. for (int i = 0; i < fNBins; i++)
  430. {
  431. for (int iCent = 0; iCent < fNCent; iCent++)
  432. {
  433. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  434. if (binCont == 0.)
  435. continue;
  436. p_Sin2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  437. }
  438. }
  439. delete prof;
  440. // WEST
  441. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  442. for (int i = 0; i < fNBins; i++)
  443. {
  444. for (int iCent = 0; iCent < fNCent; iCent++)
  445. {
  446. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  447. if (binCont == 0.)
  448. continue;
  449. p_Cos2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  450. }
  451. }
  452. delete prof;
  453. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  454. for (int i = 0; i < fNBins; i++)
  455. {
  456. for (int iCent = 0; iCent < fNCent; iCent++)
  457. {
  458. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  459. if (binCont == 0.)
  460. continue;
  461. p_Sin2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  462. }
  463. }
  464. delete prof;
  465. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  466. for (int i = 0; i < fNBins; i++)
  467. {
  468. for (int iCent = 0; iCent < fNCent; iCent++)
  469. {
  470. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  471. if (binCont == 0.)
  472. continue;
  473. p_Cos2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  474. }
  475. }
  476. delete prof;
  477. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
  478. for (int i = 0; i < fNBins; i++)
  479. {
  480. for (int iCent = 0; iCent < fNCent; iCent++)
  481. {
  482. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  483. if (binCont == 0.)
  484. continue;
  485. p_Sin2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  486. }
  487. }
  488. delete prof;
  489. // for 3rd harmonics
  490. // EAST
  491. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  492. for (int i = 0; i < fNBins; i++)
  493. {
  494. for (int iCent = 0; iCent < fNCent; iCent++)
  495. {
  496. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  497. if (binCont == 0.)
  498. continue;
  499. p_Cos3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  500. }
  501. }
  502. delete prof;
  503. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  504. for (int i = 0; i < fNBins; i++)
  505. {
  506. for (int iCent = 0; iCent < fNCent; iCent++)
  507. {
  508. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  509. if (binCont == 0.)
  510. continue;
  511. p_Sin3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  512. }
  513. }
  514. delete prof;
  515. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  516. for (int i = 0; i < fNBins; i++)
  517. {
  518. for (int iCent = 0; iCent < fNCent; iCent++)
  519. {
  520. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  521. if (binCont == 0.)
  522. continue;
  523. p_Cos3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  524. }
  525. }
  526. delete prof;
  527. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  528. for (int i = 0; i < fNBins; i++)
  529. {
  530. for (int iCent = 0; iCent < fNCent; iCent++)
  531. {
  532. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  533. if (binCont == 0.)
  534. continue;
  535. p_Sin3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  536. }
  537. }
  538. delete prof;
  539. // WEST
  540. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  541. for (int i = 0; i < fNBins; i++)
  542. {
  543. for (int iCent = 0; iCent < fNCent; iCent++)
  544. {
  545. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  546. if (binCont == 0.)
  547. continue;
  548. p_Cos3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  549. }
  550. }
  551. delete prof;
  552. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  553. for (int i = 0; i < fNBins; i++)
  554. {
  555. for (int iCent = 0; iCent < fNCent; iCent++)
  556. {
  557. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  558. if (binCont == 0.)
  559. continue;
  560. p_Sin3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  561. }
  562. }
  563. delete prof;
  564. prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  565. for (int i = 0; i < fNBins; i++)
  566. {
  567. for (int iCent = 0; iCent < fNCent; iCent++)
  568. {
  569. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  570. if (binCont == 0.)
  571. continue;
  572. p_Cos3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  573. }
  574. }
  575. delete prof;
  576. prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
  577. for (int i = 0; i < fNBins; i++)
  578. {
  579. for (int iCent = 0; iCent < fNCent; iCent++)
  580. {
  581. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  582. if (binCont == 0.)
  583. continue;
  584. p_Sin3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
  585. }
  586. }
  587. delete prof;
  588. }
  589. }
  590. }
  591. //Manage gain correction
  592. TProfile2D *BBCgain_East[fArraySize];
  593. TProfile2D *BBCgain_West[fArraySize];
  594. TFile *fiGain = new TFile(gainBBC, "read");
  595. fiGain->cd();
  596. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  597. {
  598. BBCgain_East[iVtx] = (TProfile2D *)fiGain->Get(Form("p_gainBBC_East_Vtx%i", iVtx));
  599. BBCgain_West[iVtx] = (TProfile2D *)fiGain->Get(Form("p_gainBBC_West_Vtx%i", iVtx));
  600. }
  601. //Manage recentering BBC
  602. std::map<Double_t, Double_t> p_BBC_Qx_Full_EP[fNharmonics][fArraySize][fNCent];
  603. std::map<Double_t, Double_t> p_BBC_Qy_Full_EP[fNharmonics][fArraySize][fNCent];
  604. std::map<Double_t, Double_t> p_BBC_Qx_East_EP[fNharmonics][fArraySize][fNCent];
  605. std::map<Double_t, Double_t> p_BBC_Qy_East_EP[fNharmonics][fArraySize][fNCent];
  606. std::map<Double_t, Double_t> p_BBC_Qx_West_EP[fNharmonics][fArraySize][fNCent];
  607. std::map<Double_t, Double_t> p_BBC_Qy_West_EP[fNharmonics][fArraySize][fNCent];
  608. TFile *fiReCentBBC = new TFile(reCentFileBBC, "read");
  609. std::cout << "Reading file for BBC recentering..." << std::endl;
  610. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  611. {
  612. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  613. {
  614. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
  615. prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
  616. for (int i = 0; i < fNBins; i++)
  617. {
  618. for (int iCent = 0; iCent < fNCent; iCent++)
  619. {
  620. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  621. if (binCont == 0.)
  622. continue;
  623. p_BBC_Qx_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  624. }
  625. }
  626. delete prof;
  627. prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
  628. for (int i = 0; i < fNBins; i++)
  629. {
  630. for (int iCent = 0; iCent < fNCent; iCent++)
  631. {
  632. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  633. if (binCont == 0.)
  634. continue;
  635. p_BBC_Qy_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  636. }
  637. }
  638. delete prof;
  639. prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
  640. for (int i = 0; i < fNBins; i++)
  641. {
  642. for (int iCent = 0; iCent < fNCent; iCent++)
  643. {
  644. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  645. if (binCont == 0.)
  646. continue;
  647. p_BBC_Qx_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  648. }
  649. }
  650. delete prof;
  651. prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
  652. for (int i = 0; i < fNBins; i++)
  653. {
  654. for (int iCent = 0; iCent < fNCent; iCent++)
  655. {
  656. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  657. if (binCont == 0.)
  658. continue;
  659. p_BBC_Qy_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  660. }
  661. }
  662. delete prof;
  663. prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
  664. for (int i = 0; i < fNBins; i++)
  665. {
  666. for (int iCent = 0; iCent < fNCent; iCent++)
  667. {
  668. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  669. if (binCont == 0.)
  670. continue;
  671. p_BBC_Qx_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  672. }
  673. }
  674. delete prof;
  675. prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
  676. for (int i = 0; i < fNBins; i++)
  677. {
  678. for (int iCent = 0; iCent < fNCent; iCent++)
  679. {
  680. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  681. if (binCont == 0.)
  682. continue;
  683. p_BBC_Qy_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
  684. }
  685. }
  686. delete prof;
  687. }
  688. }
  689. //Manage ShiftCorr BBC
  690. std::map<Double_t, Double_t> p_BBC_Cos_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  691. std::map<Double_t, Double_t> p_BBC_Sin_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  692. std::map<Double_t, Double_t> p_BBC_Cos_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  693. std::map<Double_t, Double_t> p_BBC_Sin_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  694. std::map<Double_t, Double_t> p_BBC_Cos_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  695. std::map<Double_t, Double_t> p_BBC_Sin_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
  696. TFile *fiShiftBBC = new TFile(shiftFileBBC, "read");
  697. fiShiftBBC->cd();
  698. std::cout << "Reading file for BBC shift corrections..." << std::endl;
  699. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  700. {
  701. for (int iVtx = 0; iVtx < fArraySize; iVtx++)
  702. {
  703. for (int iShift = 0; iShift < NShiftOrderMax * 3; iShift++)
  704. {
  705. std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
  706. prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_Full_EP%i_vtx%i", iShift, 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_Cos_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  715. }
  716. }
  717. delete prof;
  718. prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_Full_EP%i_vtx%i", iShift, 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_Sin_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  727. }
  728. }
  729. delete prof;
  730. prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_East_EP%i_vtx%i", iShift, 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_Cos_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  739. }
  740. }
  741. delete prof;
  742. prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
  743. for (int i = 0; i < fNBins; i++)
  744. {
  745. for (int iCent = 0; iCent < fNCent; iCent++)
  746. {
  747. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  748. if (binCont == 0.)
  749. continue;
  750. p_BBC_Sin_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  751. }
  752. }
  753. delete prof;
  754. prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
  755. for (int i = 0; i < fNBins; i++)
  756. {
  757. for (int iCent = 0; iCent < fNCent; iCent++)
  758. {
  759. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  760. if (binCont == 0.)
  761. continue;
  762. p_BBC_Cos_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  763. }
  764. }
  765. delete prof;
  766. prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
  767. for (int i = 0; i < fNBins; i++)
  768. {
  769. for (int iCent = 0; iCent < fNCent; iCent++)
  770. {
  771. binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
  772. if (binCont == 0.)
  773. continue;
  774. p_BBC_Sin_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
  775. }
  776. }
  777. delete prof;
  778. }
  779. }
  780. }
  781. //Eventwise parameters
  782. //hHisto[0] - before event selection
  783. //hHisto[1] - after event selection
  784. TH2D *hVtxXY[2];
  785. TH1D *hVtxZ[2];
  786. TH1D *hVpdZ[2];
  787. TH1D *hVtxVpdZ[2];
  788. TProfile *pVtxZvsRun[2];
  789. TProfile *pVtxXvsRun[2];
  790. TProfile *pVtxYvsRun[2];
  791. TProfile *pVpdZvsRun[2];
  792. TProfile *pVtxZVpdZvsRun[2];
  793. TProfile *pRefMultvsRun[2];
  794. TProfile *pAdcvsRun[2];
  795. TProfile *pTofMultvsRun[2];
  796. TProfile *pTofMatchvsRun[2];
  797. TH2D *hNevVsRun = new TH2D(Form("hNevVsRun"),Form("hNevVsRun"),fNBins, runId_min - 0.5, runId_max + 0.5,16,-0.5,15.5);
  798. //BBC EP related parameters
  799. // Raw
  800. TH2D* hBBCFullQxRaw[fNharmonics+1];
  801. TH2D* hBBCFullQyRaw[fNharmonics+1];
  802. TH2D* hBBCEastQxRaw[fNharmonics+1];
  803. TH2D* hBBCEastQyRaw[fNharmonics+1];
  804. TH2D* hBBCWestQxRaw[fNharmonics+1];
  805. TH2D* hBBCWestQyRaw[fNharmonics+1];
  806. TH2D* hBBCFullPsiRaw[fNharmonics+1];
  807. TH2D* hBBCEastPsiRaw[fNharmonics+1];
  808. TH2D* hBBCWestPsiRaw[fNharmonics+1];
  809. // w/ Gain correction
  810. TH2D* hBBCFullQxGain[fNharmonics+1];
  811. TH2D* hBBCFullQyGain[fNharmonics+1];
  812. TH2D* hBBCEastQxGain[fNharmonics+1];
  813. TH2D* hBBCEastQyGain[fNharmonics+1];
  814. TH2D* hBBCWestQxGain[fNharmonics+1];
  815. TH2D* hBBCWestQyGain[fNharmonics+1];
  816. TH2D* hBBCFullPsiGain[fNharmonics+1];
  817. TH2D* hBBCEastPsiGain[fNharmonics+1];
  818. TH2D* hBBCWestPsiGain[fNharmonics+1];
  819. // w/ Gain + recentering
  820. TH2D* hBBCFullQxRec[fNharmonics+1];
  821. TH2D* hBBCFullQyRec[fNharmonics+1];
  822. TH2D* hBBCEastQxRec[fNharmonics+1];
  823. TH2D* hBBCEastQyRec[fNharmonics+1];
  824. TH2D* hBBCWestQxRec[fNharmonics+1];
  825. TH2D* hBBCWestQyRec[fNharmonics+1];
  826. TH2D* hBBCFullPsiRec[fNharmonics+1];
  827. TH2D* hBBCEastPsiRec[fNharmonics+1];
  828. TH2D* hBBCWestPsiRec[fNharmonics+1];
  829. // w/ Gain + recentering + shift
  830. TH2D* hBBCFullPsiShift[fNharmonics+1];
  831. TH2D* hBBCEastPsiShift[fNharmonics+1];
  832. TH2D* hBBCWestPsiShift[fNharmonics+1];
  833. //ZDC EP related parameters
  834. // Raw
  835. TH2D* hZDCFullQxRaw;
  836. TH2D* hZDCFullQyRaw;
  837. TH2D* hZDCEastQxRaw;
  838. TH2D* hZDCEastQyRaw;
  839. TH2D* hZDCWestQxRaw;
  840. TH2D* hZDCWestQyRaw;
  841. TH2D* hZDCFullPsiRaw;
  842. TH2D* hZDCEastPsiRaw;
  843. TH2D* hZDCWestPsiRaw;
  844. // w/ recentering
  845. TH2D* hZDCFullQxRec;
  846. TH2D* hZDCFullQyRec;
  847. TH2D* hZDCEastQxRec;
  848. TH2D* hZDCEastQyRec;
  849. TH2D* hZDCWestQxRec;
  850. TH2D* hZDCWestQyRec;
  851. TH2D* hZDCFullPsiRec;
  852. TH2D* hZDCEastPsiRec;
  853. TH2D* hZDCWestPsiRec;
  854. // w/ recentering + shift
  855. TH2D* hZDCFullPsiShift;
  856. TH2D* hZDCEastPsiShift;
  857. TH2D* hZDCWestPsiShift;
  858. //TPC EP related parameters
  859. // Raw
  860. TH2D* hTPCFullQxRaw[fNharmonics];
  861. TH2D* hTPCFullQyRaw[fNharmonics];
  862. TH2D* hTPCEastQxRaw[fNharmonics];
  863. TH2D* hTPCEastQyRaw[fNharmonics];
  864. TH2D* hTPCWestQxRaw[fNharmonics];
  865. TH2D* hTPCWestQyRaw[fNharmonics];
  866. TH2D* hTPCFullPsiRaw[fNharmonics];
  867. TH2D* hTPCEastPsiRaw[fNharmonics];
  868. TH2D* hTPCWestPsiRaw[fNharmonics];
  869. // w/ recentering
  870. TH2D* hTPCFullQxRec[fNharmonics];
  871. TH2D* hTPCFullQyRec[fNharmonics];
  872. TH2D* hTPCEastQxRec[fNharmonics];
  873. TH2D* hTPCEastQyRec[fNharmonics];
  874. TH2D* hTPCWestQxRec[fNharmonics];
  875. TH2D* hTPCWestQyRec[fNharmonics];
  876. TH2D* hTPCFullPsiRec[fNharmonics];
  877. TH2D* hTPCEastPsiRec[fNharmonics];
  878. TH2D* hTPCWestPsiRec[fNharmonics];
  879. // w/ recentering + shift
  880. TH2D* hTPCFullPsiShift[fNharmonics];
  881. TH2D* hTPCEastPsiShift[fNharmonics];
  882. TH2D* hTPCWestPsiShift[fNharmonics];
  883. TProfile *pQx2TPCvsRun[2]; // 0 - East, 1 - West
  884. TProfile *pQy2TPCvsRun[2]; // 0 - East, 1 - West
  885. TProfile *pQx3TPCvsRun[2]; // 0 - East, 1 - West
  886. TProfile *pQy3TPCvsRun[2]; // 0 - East, 1 - West
  887. TProfile *pQx2BBCvsRun[2]; // 0 - East, 1 - West
  888. TProfile *pQy2BBCvsRun[2]; // 0 - East, 1 - West
  889. TProfile *pQx1ZDCvsRun[2]; // 0 - East, 1 - West
  890. TProfile *pQy1ZDCvsRun[2]; // 0 - East, 1 - West
  891. TProfile *pCos2TPCvsRun[2*fNharmonics][2];
  892. TProfile *pSin2TPCvsRun[2*fNharmonics][2];
  893. TProfile *pCos3TPCvsRun[2*fNharmonics][2];
  894. TProfile *pSin3TPCvsRun[2*fNharmonics][2];
  895. TProfile *pCosBBCvsRun[2*fNharmonics][2];
  896. TProfile *pSinBBCvsRun[2*fNharmonics][2];
  897. TProfile *pCosZDCvsRun[2];
  898. TProfile *pSinZDCvsRun[2];
  899. //Trackwise parameters
  900. //hHisto[0] - before track selection (but after event selection)
  901. //hHisto[1] - after track selection
  902. TH1D *hNhitsFit[2];
  903. TH1D *hNhitsPoss[2];
  904. TH1D *hNhitsRatio[2];
  905. TH1D *hPt[2];
  906. TH1D *hPtot[2];
  907. TH1D *hPhi[2];
  908. TH1D *hEta[2];
  909. TH1D *hDCA[2];
  910. TProfile *pPtvsRun[2];
  911. TProfile *pEtavsRun[2];
  912. TProfile *pPhivsRun[2];
  913. TProfile *pNhitsFitvsRun[2];
  914. TProfile *pNhitsPossvsRun[2];
  915. TProfile *pDCAvsRun[2];
  916. TProfile *pdEdxvsRun[2];
  917. TProfile *pNtrPrimvsRun[2];
  918. TProfile *pNtrGlobvsRun[2];
  919. //PID related parameters
  920. //hHisto[0] - before any PID (but after track selection)
  921. //hHisto[1] - after PID
  922. //second [i] is for particle species: 0-all, 1-pion, 2-kaon, 3-proton
  923. TH2D *hnSigM2_lowpt[2];
  924. TH2D *hnSigM2_highpt[2];
  925. TH2D *hdEdxQp[2][4];
  926. TH2D *hM2Qp[2][4];
  927. TH2D *hNsigQp[2][4];
  928. TH2D *hM2pt[2][4];
  929. TH2D *hNsigpt[2][4];
  930. TProfile2D *BBCgain[2][24];
  931. TH2D *h_pt_spectra[4][16]; //for <pT> calc 0 - pions, 1 - kaons, 2 - protons, 3 - charged hadrons
  932. //Initialization
  933. //ZDC
  934. // Raw
  935. hZDCFullQxRaw = new TH2D(Form("hZDCFullQxRaw%i",0), Form("hZDCFullQxRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  936. hZDCFullQyRaw = new TH2D(Form("hZDCFullQyRaw%i",0), Form("hZDCFullQyRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  937. hZDCEastQxRaw = new TH2D(Form("hZDCEastQxRaw%i",0), Form("hZDCEastQxRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  938. hZDCEastQyRaw = new TH2D(Form("hZDCEastQyRaw%i",0), Form("hZDCEastQyRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  939. hZDCWestQxRaw = new TH2D(Form("hZDCWestQxRaw%i",0), Form("hZDCWestQxRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  940. hZDCWestQyRaw = new TH2D(Form("hZDCWestQyRaw%i",0), Form("hZDCWestQyRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  941. hZDCFullPsiRaw = new TH2D(Form("hZDCFullPsiRaw%i",0), Form("hZDCFullPsiRaw%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  942. hZDCEastPsiRaw = new TH2D(Form("hZDCEastPsiRaw%i",0), Form("hZDCEastPsiRaw%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  943. hZDCWestPsiRaw = new TH2D(Form("hZDCWestPsiRaw%i",0), Form("hZDCWestPsiRaw%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  944. // Rec
  945. hZDCFullQxRec = new TH2D(Form("hZDCFullQxRec%i",0), Form("hZDCFullQxRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  946. hZDCFullQyRec = new TH2D(Form("hZDCFullQyRec%i",0), Form("hZDCFullQyRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  947. hZDCEastQxRec = new TH2D(Form("hZDCEastQxRec%i",0), Form("hZDCEastQxRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  948. hZDCEastQyRec = new TH2D(Form("hZDCEastQyRec%i",0), Form("hZDCEastQyRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  949. hZDCWestQxRec = new TH2D(Form("hZDCWestQxRec%i",0), Form("hZDCWestQxRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  950. hZDCWestQyRec = new TH2D(Form("hZDCWestQyRec%i",0), Form("hZDCWestQyRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
  951. hZDCFullPsiRec = new TH2D(Form("hZDCFullPsiRec%i",0), Form("hZDCFullPsiRec%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  952. hZDCEastPsiRec = new TH2D(Form("hZDCEastPsiRec%i",0), Form("hZDCEastPsiRec%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  953. hZDCWestPsiRec = new TH2D(Form("hZDCWestPsiRec%i",0), Form("hZDCWestPsiRec%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  954. // Shift
  955. hZDCFullPsiShift = new TH2D(Form("hZDCFullPsiShift%i",0), Form("hZDCFullPsiShift%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  956. hZDCEastPsiShift = new TH2D(Form("hZDCEastPsiShift%i",0), Form("hZDCEastPsiShift%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  957. hZDCWestPsiShift = new TH2D(Form("hZDCWestPsiShift%i",0), Form("hZDCWestPsiShift%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
  958. //BBC
  959. for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
  960. {
  961. // Raw
  962. hBBCFullQxRaw[iHarm] = new TH2D(Form("hBBCFullQxRaw%i",iHarm), Form("hBBCFullQxRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  963. hBBCFullQyRaw[iHarm] = new TH2D(Form("hBBCFullQyRaw%i",iHarm), Form("hBBCFullQyRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  964. hBBCEastQxRaw[iHarm] = new TH2D(Form("hBBCEastQxRaw%i",iHarm), Form("hBBCEastQxRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  965. hBBCEastQyRaw[iHarm] = new TH2D(Form("hBBCEastQyRaw%i",iHarm), Form("hBBCEastQyRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  966. hBBCWestQxRaw[iHarm] = new TH2D(Form("hBBCWestQxRaw%i",iHarm), Form("hBBCWestQxRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  967. hBBCWestQyRaw[iHarm] = new TH2D(Form("hBBCWestQyRaw%i",iHarm), Form("hBBCWestQyRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  968. hBBCFullPsiRaw[iHarm] = new TH2D(Form("hBBCFullPsiRaw%i",iHarm), Form("hBBCFullPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  969. hBBCEastPsiRaw[iHarm] = new TH2D(Form("hBBCEastPsiRaw%i",iHarm), Form("hBBCEastPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  970. hBBCWestPsiRaw[iHarm] = new TH2D(Form("hBBCWestPsiRaw%i",iHarm), Form("hBBCWestPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  971. // Gain
  972. hBBCFullQxGain[iHarm] = new TH2D(Form("hBBCFullQxGain%i",iHarm), Form("hBBCFullQxGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  973. hBBCFullQyGain[iHarm] = new TH2D(Form("hBBCFullQyGain%i",iHarm), Form("hBBCFullQyGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  974. hBBCEastQxGain[iHarm] = new TH2D(Form("hBBCEastQxGain%i",iHarm), Form("hBBCEastQxGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  975. hBBCEastQyGain[iHarm] = new TH2D(Form("hBBCEastQyGain%i",iHarm), Form("hBBCEastQyGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  976. hBBCWestQxGain[iHarm] = new TH2D(Form("hBBCWestQxGain%i",iHarm), Form("hBBCWestQxGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  977. hBBCWestQyGain[iHarm] = new TH2D(Form("hBBCWestQyGain%i",iHarm), Form("hBBCWestQyGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  978. hBBCFullPsiGain[iHarm] = new TH2D(Form("hBBCFullPsiGain%i",iHarm), Form("hBBCFullPsiGain%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  979. hBBCEastPsiGain[iHarm] = new TH2D(Form("hBBCEastPsiGain%i",iHarm), Form("hBBCEastPsiGain%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  980. hBBCWestPsiGain[iHarm] = new TH2D(Form("hBBCWestPsiGain%i",iHarm), Form("hBBCWestPsiGain%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  981. // Rec
  982. hBBCFullQxRec[iHarm] = new TH2D(Form("hBBCFullQxRec%i",iHarm), Form("hBBCFullQxRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  983. hBBCFullQyRec[iHarm] = new TH2D(Form("hBBCFullQyRec%i",iHarm), Form("hBBCFullQyRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  984. hBBCEastQxRec[iHarm] = new TH2D(Form("hBBCEastQxRec%i",iHarm), Form("hBBCEastQxRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  985. hBBCEastQyRec[iHarm] = new TH2D(Form("hBBCEastQyRec%i",iHarm), Form("hBBCEastQyRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  986. hBBCWestQxRec[iHarm] = new TH2D(Form("hBBCWestQxRec%i",iHarm), Form("hBBCWestQxRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  987. hBBCWestQyRec[iHarm] = new TH2D(Form("hBBCWestQyRec%i",iHarm), Form("hBBCWestQyRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
  988. hBBCFullPsiRec[iHarm] = new TH2D(Form("hBBCFullPsiRec%i",iHarm), Form("hBBCFullPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  989. hBBCEastPsiRec[iHarm] = new TH2D(Form("hBBCEastPsiRec%i",iHarm), Form("hBBCEastPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  990. hBBCWestPsiRec[iHarm] = new TH2D(Form("hBBCWestPsiRec%i",iHarm), Form("hBBCWestPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  991. // Shift
  992. hBBCFullPsiShift[iHarm] = new TH2D(Form("hBBCFullPsiShift%i",iHarm), Form("hBBCFullPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  993. hBBCEastPsiShift[iHarm] = new TH2D(Form("hBBCEastPsiShift%i",iHarm), Form("hBBCEastPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  994. hBBCWestPsiShift[iHarm] = new TH2D(Form("hBBCWestPsiShift%i",iHarm), Form("hBBCWestPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  995. }
  996. //TPC
  997. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  998. {
  999. // Raw
  1000. hTPCFullQxRaw[iHarm] = new TH2D(Form("hTPCFullQxRaw%i",iHarm), Form("hTPCFullQxRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1001. hTPCFullQyRaw[iHarm] = new TH2D(Form("hTPCFullQyRaw%i",iHarm), Form("hTPCFullQyRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1002. hTPCEastQxRaw[iHarm] = new TH2D(Form("hTPCEastQxRaw%i",iHarm), Form("hTPCEastQxRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1003. hTPCEastQyRaw[iHarm] = new TH2D(Form("hTPCEastQyRaw%i",iHarm), Form("hTPCEastQyRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1004. hTPCWestQxRaw[iHarm] = new TH2D(Form("hTPCWestQxRaw%i",iHarm), Form("hTPCWestQxRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1005. hTPCWestQyRaw[iHarm] = new TH2D(Form("hTPCWestQyRaw%i",iHarm), Form("hTPCWestQyRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1006. hTPCFullPsiRaw[iHarm] = new TH2D(Form("hTPCFullPsiRaw%i",iHarm), Form("hTPCFullPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1007. hTPCEastPsiRaw[iHarm] = new TH2D(Form("hTPCEastPsiRaw%i",iHarm), Form("hTPCEastPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1008. hTPCWestPsiRaw[iHarm] = new TH2D(Form("hTPCWestPsiRaw%i",iHarm), Form("hTPCWestPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1009. // Rec
  1010. hTPCFullQxRec[iHarm] = new TH2D(Form("hTPCFullQxRec%i",iHarm), Form("hTPCFullQxRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1011. hTPCFullQyRec[iHarm] = new TH2D(Form("hTPCFullQyRec%i",iHarm), Form("hTPCFullQyRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1012. hTPCEastQxRec[iHarm] = new TH2D(Form("hTPCEastQxRec%i",iHarm), Form("hTPCEastQxRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1013. hTPCEastQyRec[iHarm] = new TH2D(Form("hTPCEastQyRec%i",iHarm), Form("hTPCEastQyRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1014. hTPCWestQxRec[iHarm] = new TH2D(Form("hTPCWestQxRec%i",iHarm), Form("hTPCWestQxRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1015. hTPCWestQyRec[iHarm] = new TH2D(Form("hTPCWestQyRec%i",iHarm), Form("hTPCWestQyRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
  1016. hTPCFullPsiRec[iHarm] = new TH2D(Form("hTPCFullPsiRec%i",iHarm), Form("hTPCFullPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1017. hTPCEastPsiRec[iHarm] = new TH2D(Form("hTPCEastPsiRec%i",iHarm), Form("hTPCEastPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1018. hTPCWestPsiRec[iHarm] = new TH2D(Form("hTPCWestPsiRec%i",iHarm), Form("hTPCWestPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1019. // Shift
  1020. hTPCFullPsiShift[iHarm] = new TH2D(Form("hTPCFullPsiShift%i",iHarm), Form("hTPCFullPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1021. hTPCEastPsiShift[iHarm] = new TH2D(Form("hTPCEastPsiShift%i",iHarm), Form("hTPCEastPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1022. hTPCWestPsiShift[iHarm] = new TH2D(Form("hTPCWestPsiShift%i",iHarm), Form("hTPCWestPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
  1023. }
  1024. for (int i = 0; i < 2; i++)
  1025. {
  1026. //Eventwise
  1027. hVtxXY[i] = new TH2D(Form("hVtxXY%i", i), Form("vertex XY %i;Vtx_{X}, [cm];Vtx_{Y}, [cm]", i), 200, -5., 5., 200, -5., 5.);
  1028. hVtxZ[i] = new TH1D(Form("hVtxZ%i", i), Form("vertex Z %i;Vtx_{Z}, [cm];N_{counts}", i), 200, -40., 40.);
  1029. hVpdZ[i] = new TH1D(Form("hVpdZ%i", i), Form("vertex Z %i;Vpd_{Z}, [cm];N_{counts}", i), 200, -40., 40.);
  1030. hVtxVpdZ[i] = new TH1D(Form("hVtxVpdZ%i", i), Form("VtxZ - VpdVz %i;Vtx_{Z}-vpdV_{Z}, [cm];N_{counts}", i), 200, -5., 5.);
  1031. pVtxZvsRun[i] = new TProfile(Form("pVtxZvsRun%i",i),Form("pVtxZvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1032. pVtxXvsRun[i] = new TProfile(Form("pVtxXvsRun%i",i),Form("pVtxXvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1033. pVtxYvsRun[i] = new TProfile(Form("pVtxYvsRun%i",i),Form("pVtxYvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1034. pVpdZvsRun[i] = new TProfile(Form("pVpdZvsRun%i",i),Form("pVpdZvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1035. pVtxZVpdZvsRun[i] = new TProfile(Form("pVtxZVpdZvsRun%i",i),Form("pVtxZVpdZvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1036. pRefMultvsRun[i] = new TProfile(Form("pRefMultvsRun%i",i),Form("pRefMultvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1037. pAdcvsRun[i] = new TProfile(Form("pAdcvsRun%i",i),Form("pAdcvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1038. pTofMultvsRun[i] = new TProfile(Form("pTofMultvsRun%i",i),Form("pTofMultvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1039. pTofMatchvsRun[i] = new TProfile(Form("pTofMatchvsRun%i",i),Form("pTofMatchvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1040. //Trackwise
  1041. hNhitsFit[i] = new TH1D(Form("hNhitsFit%i", i), Form("N_{hist}^{Fit} %i;N_{hist}^{Fit};N_{counts}", i), 50, 0., 50.);
  1042. hNhitsPoss[i] = new TH1D(Form("hNhitsPoss%i", i), Form("N_{hist}^{Poss} %i;N_{hist}^{Poss};N_{counts}", i), 50, 0., 50.);
  1043. hNhitsRatio[i] = new TH1D(Form("hNhitsRatio%i", i), Form("N_{hist}^{Fit}/N_{hist}^{Poss} %i;N_{hist}^{Fit}/N_{hist}^{Poss};N_{counts}", i), 50, 0., 1.);
  1044. hPt[i] = new TH1D(Form("hPt%i", i), Form("p_{T} %i;p_{T}, [GeV/c];N_{counts}", i), 600, 0., 6.);
  1045. hPtot[i] = new TH1D(Form("hPtot%i", i), Form("p_{tot} %i;p, [GeV/c];N_{counts}", i), 600, 0., 12.);
  1046. hPhi[i] = new TH1D(Form("hPhi%i", i), Form("#varphi %i;#varphi, [rad];N_{counts}", i), 600, 0., 6.);
  1047. hEta[i] = new TH1D(Form("hEta%i", i), Form("#eta %i;#eta;N_{counts}", i), 400, -2., 2.);
  1048. hDCA[i] = new TH1D(Form("hDCA%i", i), Form("|DCA| %i;dca, [cm];N_{counts}", i), 200, 0., 5.);
  1049. pPtvsRun[i] = new TProfile(Form("pPtvsRun%i",i),Form("pPtvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1050. pEtavsRun[i] = new TProfile(Form("pEtavsRun%i",i),Form("pEtavsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1051. pPhivsRun[i] = new TProfile(Form("pPhivsRun%i",i),Form("pPhivsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1052. pNhitsFitvsRun[i] = new TProfile(Form("pNhitsFitvsRun%i",i),Form("pNhitsFitvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1053. pNhitsPossvsRun[i] = new TProfile(Form("pNhitsPossvsRun%i",i),Form("pNhitsPossvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1054. pDCAvsRun[i] = new TProfile(Form("pDCAvsRun%i",i),Form("pDCAvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1055. pdEdxvsRun[i] = new TProfile(Form("pdEdxvsRun%i",i),Form("pdEdxvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1056. pNtrPrimvsRun[i] = new TProfile(Form("pNtrPrimvsRun%i",i),Form("pNtrPrimvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1057. pNtrGlobvsRun[i] = new TProfile(Form("pNtrGlobvsRun%i",i),Form("pNtrGlobvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1058. pQx2TPCvsRun[i] = new TProfile(Form("pQx2TPCvsRun%i",i),Form("pQx2TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1059. pQy2TPCvsRun[i] = new TProfile(Form("pQy2TPCvsRun%i",i),Form("pQy2TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1060. pQx3TPCvsRun[i] = new TProfile(Form("pQx3TPCvsRun%i",i),Form("pQx3TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1061. pQy3TPCvsRun[i] = new TProfile(Form("pQy3TPCvsRun%i",i),Form("pQy3TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1062. pQx2BBCvsRun[i] = new TProfile(Form("pQx2BBCvsRun%i",i),Form("pQx2BBCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1063. pQy2BBCvsRun[i] = new TProfile(Form("pQy2BBCvsRun%i",i),Form("pQy2BBCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1064. for (int iHarm=0; iHarm<2*fNharmonics; iHarm++)
  1065. {
  1066. pCos2TPCvsRun[iHarm][i] = new TProfile(Form("pCos2_%iTPCvsRun%i",iHarm,i),Form("pCos2_%iTPCvsRun%i",iHarm,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1067. pSin2TPCvsRun[iHarm][i] = new TProfile(Form("pSin2_%iTPCvsRun%i",iHarm,i),Form("pSin2_%iTPCvsRun%i",iHarm,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1068. pCos3TPCvsRun[iHarm][i] = new TProfile(Form("pCos3_%iTPCvsRun%i",iHarm,i),Form("pCos3_%iTPCvsRun%i",iHarm,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1069. pSin3TPCvsRun[iHarm][i] = new TProfile(Form("pSin3_%iTPCvsRun%i",iHarm,i),Form("pSin3_%iTPCvsRun%i",iHarm,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1070. pCosBBCvsRun[iHarm][i] = new TProfile(Form("pCos%iBBCvsRun%i",iHarm,i),Form("pCos%iBBCvsRun%i",iHarm,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1071. pSinBBCvsRun[iHarm][i] = new TProfile(Form("pSin%iBBCvsRun%i",iHarm,i),Form("pSin%iBBCvsRun%i",iHarm,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1072. }
  1073. pCosZDCvsRun[i] = new TProfile(Form("pCos%iZDCvsRun%i",0,i),Form("pCos%iZDCvsRun%i",0,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1074. pSinZDCvsRun[i] = new TProfile(Form("pSin%iZDCvsRun%i",0,i),Form("pSin%iZDCvsRun%i",0,i),fNBins, runId_min - 0.5, runId_max + 0.5);
  1075. //PID
  1076. hnSigM2_lowpt[i] = new TH2D(Form("hnSigM2_lowpt%i", i), Form("n#sigma_{#pi} vs M^{2} for 0.5 < p_{T} < 0.7 GeV/c %i;n#sigma_{#pi};m^{2} [GeV/c^{2}]^{2}", i), 600, -16., 16., 600, -0.8, 1.8);
  1077. hnSigM2_highpt[i] = new TH2D(Form("hnSigM2_highpt%i", i), Form("n#sigma_{#pi} vs M^{2} for 2.2 < p_{T} < 2.4 GeV/c %i;n#sigma_{#pi};m^{2} [GeV/c^{2}]^{2}", i), 650, -7., 6., 600, -0.8, 1.8);
  1078. for (int iPID = 0; iPID < 4; iPID++)
  1079. {
  1080. hdEdxQp[i][iPID] = new TH2D(Form("hdEdxQp%i%i", i, iPID), Form("dEdx vs Q*p_{tot} %i %i;Q*p, [GeV/c];dEdx, [a.u.]", i, iPID), 1000, -5., 5., 1000, 0., 2e-5);
  1081. hM2Qp[i][iPID] = new TH2D(Form("hM2Qp%i%i", i, iPID), Form("m^{2} vs Q*p_{tot} %i %i;Q*p, [GeV/c];m^{2}, [GeV/c^{2}]^{2}", i, iPID), 1000, -5., 5., 1000, -0.8, 1.8);
  1082. hNsigQp[i][iPID] = new TH2D(Form("hNsigQp%i%i", i, iPID), Form("n#sigma_{#pi} vs Q*p_{tot} %i %i;Q*p, [GeV/c];n#sigma_{#pi}", i, iPID), 1000, -5., 5., 600, -3., 3.);
  1083. hM2pt[i][iPID] = new TH2D(Form("hM2pt%i%i", i, iPID), Form("m^{2} vs p_{T} %i %i;p_{T}, [GeV/c];m^{2}, [GeV/c^{2}]^{2}", i, iPID), 1000, 0., 5., 1000, -0.8, 1.8);
  1084. hNsigpt[i][iPID] = new TH2D(Form("hNsigpt%i%i", i, iPID), Form("n#sigma_{#pi} vs p_{T} %i %i;p_{T}, [GeV/c];n#sigma_{#pi}", i, iPID), 1000, 0., 5., 600, -3., 3.);
  1085. }
  1086. }
  1087. for (int iPID = 0; iPID < 4; iPID++)
  1088. {
  1089. for (int iCent=0;iCent<16;iCent++)
  1090. {
  1091. h_pt_spectra[iPID][iCent] = new TH2D(Form("h_pt_spectra_particle%i_cent%i",iPID,iCent),Form("h_pt_spectra_particle%i_cent%i",iPID,iCent),fNBins, runId_min -0.5, runId_max +0.5,550,0,5.5);
  1092. // hist[iPID][cent]->Fill(run,pt);
  1093. }
  1094. }
  1095. //Counters
  1096. Int_t fUsedTracks = 0;
  1097. Int_t fQcountFull = 0;
  1098. Int_t fQcountFullEast = 0;
  1099. Int_t fQcountFullWest = 0;
  1100. Int_t fQcountWest[NEtaGaps];
  1101. Int_t fQcountEast[NEtaGaps];
  1102. //Q-vectors
  1103. TVector2 fQv2FullEP;
  1104. TVector2 fQv2FullEPRaw;
  1105. TVector2 fQv2FullSP;
  1106. TVector2 fQv3FullEP;
  1107. TVector2 fQv3FullEPRaw;
  1108. TVector2 fQv3FullSP;
  1109. TVector2 fQv2WestEP[NEtaGaps];
  1110. TVector2 fQv2EastEP[NEtaGaps];
  1111. TVector2 fQv2WestEPRaw[NEtaGaps];
  1112. TVector2 fQv2EastEPRaw[NEtaGaps];
  1113. TVector2 fQv2WestSP[NEtaGaps];
  1114. TVector2 fQv2EastSP[NEtaGaps];
  1115. TVector2 fQv3WestEP[NEtaGaps];
  1116. TVector2 fQv3EastEP[NEtaGaps];
  1117. TVector2 fQv3WestEPRaw[NEtaGaps];
  1118. TVector2 fQv3EastEPRaw[NEtaGaps];
  1119. TVector2 fQv3WestSP[NEtaGaps];
  1120. TVector2 fQv3EastSP[NEtaGaps];
  1121. TVector2 QvTMP;
  1122. TVector2 fQv2WestSP_Shift[NEtaGaps];
  1123. TVector2 fQv2EastSP_Shift[NEtaGaps];
  1124. TVector2 fQv3WestSP_Shift[NEtaGaps];
  1125. TVector2 fQv3EastSP_Shift[NEtaGaps];
  1126. TVector2 fUvSP;
  1127. Double_t weight;
  1128. Int_t VtxSign;
  1129. Double_t psiShiftedTMP;
  1130. Double_t PsiSin;
  1131. Double_t PsiCos;
  1132. Long64_t goodEventCounter = 0;
  1133. Double_t Psi2FullEP, Psi3FullEP;
  1134. Double_t Psi2EastEP[NEtaGaps], Psi3EastEP[NEtaGaps];
  1135. Double_t Psi2WestEP[NEtaGaps], Psi3WestEP[NEtaGaps];
  1136. Double_t Psi2EastSP[NEtaGaps], Psi3EastSP[NEtaGaps];
  1137. Double_t Psi2WestSP[NEtaGaps], Psi3WestSP[NEtaGaps];
  1138. Double_t Psi_ReCenter;
  1139. Int_t bin;
  1140. Double_t res2EP, res3EP, res2EPsqr, res3EPsqr;
  1141. Double_t flow2EP;
  1142. Double_t flow3EP;
  1143. Double_t res2SP, res3SP, res2SPsqr, res3SPsqr;
  1144. Double_t flow2SP;
  1145. Double_t flow3SP;
  1146. Double_t psiSP_rec;
  1147. int k = 0;
  1148. TVector2 qvBBC_East, qvBBC_West;
  1149. Int_t adcBin, shiftBin;
  1150. Double_t psiEP_rec, psiEP_BBC_Full[fNharmonics+1], psiEP_BBC_East[fNharmonics+1], psiEP_BBC_West[fNharmonics+1], shiftCos, shiftSin;
  1151. Double_t dPsi;
  1152. // Loop over events
  1153. for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
  1154. {
  1155. if ((iEvent + 1) % 1000 == 0)
  1156. {
  1157. std::cout << "Working on event #[" << (iEvent + 1)
  1158. << "/" << events2read << "]" << std::endl;
  1159. }
  1160. Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
  1161. if (!readEvent)
  1162. {
  1163. std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
  1164. break;
  1165. }
  1166. // Retrieve femtoDst
  1167. StFemtoDst *dst = femtoReader->femtoDst();
  1168. // Retrieve event information
  1169. StFemtoEvent *event = dst->event();
  1170. if (!event)
  1171. {
  1172. std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
  1173. break;
  1174. }
  1175. TVector3 pVtx = event->primaryVertex();
  1176. hVtxXY[0]->Fill(pVtx.X(), pVtx.Y());
  1177. hVtxZ[0]->Fill(pVtx.Z());
  1178. hVtxVpdZ[0]->Fill(pVtx.Z() - event->vpdVz());
  1179. hVpdZ[0]->Fill(event->vpdVz());
  1180. pVtxXvsRun[0]->Fill(event->runId(),pVtx.X());
  1181. pVtxYvsRun[0]->Fill(event->runId(),pVtx.Y());
  1182. pVtxZvsRun[0]->Fill(event->runId(),pVtx.Z());
  1183. pVpdZvsRun[0]->Fill(event->runId(),event->vpdVz());
  1184. pVtxZVpdZvsRun[0]->Fill(event->runId(),pVtx.Z() - event->vpdVz());
  1185. pRefMultvsRun[0]->Fill(event->runId(),event->refMult());
  1186. pTofMultvsRun[0]->Fill(event->runId(),event->numberOfBTofHit());
  1187. pTofMatchvsRun[0]->Fill(event->runId(),event->numberOfTofMatched());
  1188. pNtrPrimvsRun[0]->Fill(event->runId(),event->numberOfPrimaryTracks());
  1189. pNtrGlobvsRun[0]->Fill(event->runId(),event->numberOfGlobalTracks());
  1190. Int_t fAdcSum=0;
  1191. for (int i=0;i<16;i++)
  1192. {
  1193. fAdcSum += event->bbcAdcEast(i);
  1194. fAdcSum += event->bbcAdcWest(i);
  1195. }
  1196. pAdcvsRun[0]->Fill(event->runId(),fAdcSum);
  1197. // Event selection
  1198. if (!isGoodEvent(event))
  1199. continue;
  1200. VtxSign = GetVzBin(event->primaryVertex().Z());
  1201. if (VtxSign == -1)
  1202. continue;
  1203. hNevVsRun->Fill(event->runId(),event->cent16());
  1204. //Counters
  1205. fUsedTracks = 0;
  1206. fQcountFull = 0;
  1207. fQcountFullEast = 0;
  1208. fQcountFullWest = 0;
  1209. psiShiftedTMP = 0.;
  1210. PsiSin = 0.;
  1211. PsiCos = 0.;
  1212. fQv2FullEP.Set(0.,0.);
  1213. fQv2FullEPRaw.Set(0.,0.);
  1214. fQv2FullSP.Set(0.,0.);
  1215. fQv3FullEP.Set(0.,0.);
  1216. fQv3FullEPRaw.Set(0.,0.);
  1217. fQv3FullSP.Set(0.,0.);
  1218. QvTMP.Set(0.,0.);
  1219. for (int i=0; i<NEtaGaps; i++)
  1220. {
  1221. fQcountWest[i] = 0;
  1222. fQcountEast[i] = 0;
  1223. fQv2WestEP[i].Set(0.,0.);
  1224. fQv2EastEP[i].Set(0.,0.);
  1225. fQv2WestEPRaw[i].Set(0.,0.);
  1226. fQv2EastEPRaw[i].Set(0.,0.);
  1227. fQv2WestSP[i].Set(0.,0.);
  1228. fQv2EastSP[i].Set(0.,0.);
  1229. fQv3WestEP[i].Set(0.,0.);
  1230. fQv3EastEP[i].Set(0.,0.);
  1231. fQv3WestEPRaw[i].Set(0.,0.);
  1232. fQv3EastEPRaw[i].Set(0.,0.);
  1233. fQv3WestSP[i].Set(0.,0.);
  1234. fQv3EastSP[i].Set(0.,0.);
  1235. }
  1236. hVtxXY[1]->Fill(pVtx.X(), pVtx.Y());
  1237. hVtxZ[1]->Fill(pVtx.Z());
  1238. hVtxVpdZ[1]->Fill(pVtx.Z() - event->vpdVz());
  1239. hVpdZ[1]->Fill(event->vpdVz());
  1240. pVtxXvsRun[1]->Fill(event->runId(),pVtx.X());
  1241. pVtxYvsRun[1]->Fill(event->runId(),pVtx.Y());
  1242. pVtxZvsRun[1]->Fill(event->runId(),pVtx.Z());
  1243. pVpdZvsRun[1]->Fill(event->runId(),event->vpdVz());
  1244. pVtxZVpdZvsRun[1]->Fill(event->runId(),pVtx.Z() - event->vpdVz());
  1245. pRefMultvsRun[1]->Fill(event->runId(),event->refMult());
  1246. pTofMultvsRun[1]->Fill(event->runId(),event->numberOfBTofHit());
  1247. pTofMatchvsRun[1]->Fill(event->runId(),event->numberOfTofMatched());
  1248. pNtrPrimvsRun[1]->Fill(event->runId(),event->numberOfPrimaryTracks());
  1249. pNtrGlobvsRun[1]->Fill(event->runId(),event->numberOfGlobalTracks());
  1250. fAdcSum=0;
  1251. for (int i=0;i<16;i++)
  1252. {
  1253. fAdcSum += event->bbcAdcEast(i);
  1254. fAdcSum += event->bbcAdcWest(i);
  1255. }
  1256. pAdcvsRun[1]->Fill(event->runId(),fAdcSum);
  1257. if (k == 0)
  1258. {
  1259. std::cout << "East BBC tileId = 8 with phi = " << BBC_GetPhi(0, 8) << std::endl;
  1260. std::cout << "West BBC tileId = 8 with phi = " << BBC_GetPhi(1, 8) << std::endl;
  1261. k++;
  1262. }
  1263. //BBC EP
  1264. //Raw
  1265. TVector2 qRawEast[fNharmonics+1], qRawWest[fNharmonics+1], qRawFull[fNharmonics+1];
  1266. Double_t psiRawEastBBC[fNharmonics+1], psiRawWestBBC[fNharmonics+1], psiRawFullBBC[fNharmonics+1];
  1267. float nHitE_Raw[16], nHitW_Raw[16];
  1268. // 0 is for full event plane, 1 is for west only event plane, 2 is for east only event plane
  1269. //the number of hits in the bbc tile is simply the adc
  1270. for(int i=0;i<16;i++){
  1271. nHitE_Raw[i] = event->bbcAdcEast(i);
  1272. nHitW_Raw[i] = event->bbcAdcWest(i);
  1273. }
  1274. for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
  1275. {
  1276. //BBC Full
  1277. //v1*y > 0 by conventions most reasonable to the STAR event plane
  1278. //therefore the west event plane has the right sign for the event plane
  1279. //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
  1280. //The east event plane points opposite that of the west and the west finds the "correct" EP for the STAR coordinates, use West-East
  1281. Double_t QVector_RawX = 0., QVector_RawY = 0.;
  1282. Float_t eXFsum_Raw = 0., wXFsum_Raw = 0., eYFsum_Raw = 0., wYFsum_Raw = 0., eWFgt_Raw=0., wWFgt_Raw = 0.;
  1283. TVector2 QFullRaw;
  1284. for(int iTile = 0; iTile < 16; iTile++) {
  1285. eXFsum_Raw += cos((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
  1286. wXFsum_Raw += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
  1287. wYFsum_Raw += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
  1288. eYFsum_Raw += sin((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
  1289. eWFgt_Raw += nHitE_Raw[iTile];
  1290. wWFgt_Raw += nHitW_Raw[iTile];
  1291. }
  1292. QVector_RawX=(eWFgt_Raw > 0. && wWFgt_Raw > 0.) ? wXFsum_Raw/wWFgt_Raw - eXFsum_Raw/eWFgt_Raw:0.;
  1293. QVector_RawY=(eWFgt_Raw > 0. && wWFgt_Raw > 0.) ? wYFsum_Raw/wWFgt_Raw - eYFsum_Raw/eWFgt_Raw:0.;
  1294. QFullRaw.Set(QVector_RawX,QVector_RawY);
  1295. qRawFull[iHarm] = QFullRaw;
  1296. if(qRawFull[iHarm].Mod()){
  1297. psiRawFullBBC[iHarm] = qRawFull[iHarm].Phi()/(iHarm+1.);
  1298. if(psiRawFullBBC[iHarm] < 0.0) {psiRawFullBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1299. }
  1300. //BBC West
  1301. //v1*y > 0 by conventions most reasonable to the STAR event plane
  1302. //therefore the west event plane has the right sign for the event plane
  1303. //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
  1304. Double_t QVectorW_RawX = 0., QVectorW_RawY = 0.;
  1305. Float_t wXsum_Raw = 0., wYsum_Raw = 0., wWgt_Raw = 0.;
  1306. Double_t EventPlaneRawWest = 0.;
  1307. TVector2 QWestRaw;
  1308. for(int iTile = 0; iTile < 16; iTile++) {
  1309. wXsum_Raw += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
  1310. wYsum_Raw += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
  1311. wWgt_Raw += nHitW_Raw[iTile];
  1312. }
  1313. QVectorW_RawX = (wWgt_Raw > 0.0) ? wXsum_Raw/wWgt_Raw:0.0;
  1314. QVectorW_RawY = (wWgt_Raw > 0.0) ? wYsum_Raw/wWgt_Raw:0.0;
  1315. QWestRaw.Set(QVectorW_RawX,QVectorW_RawY);
  1316. qRawWest[iHarm] = QWestRaw;
  1317. if(qRawWest[iHarm].Mod()){
  1318. psiRawWestBBC[iHarm] = qRawWest[iHarm].Phi()/(iHarm+1.);
  1319. if(psiRawWestBBC[iHarm] < 0.0) {psiRawWestBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1320. }
  1321. //BBC East
  1322. //in BBC_GetPhi the fact that the tile numbering scheme is reversed about the y axis for the east BBC is accounted for
  1323. //The east event plane points opposite that of the west since the spectators are on different sides -- see west comments
  1324. Double_t QVectorE_RawX = 0., QVectorE_RawY = 0.;
  1325. Float_t eXsum_Raw = 0., eYsum_Raw = 0., eWgt_Raw = 0.;
  1326. TVector2 QEastRaw;
  1327. for(int iTile = 0; iTile < 16; iTile++) {
  1328. eXsum_Raw += cos(BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
  1329. eYsum_Raw += sin(BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
  1330. eWgt_Raw += nHitE_Raw[iTile];
  1331. }
  1332. QVectorE_RawX = (eWgt_Raw > 0.0) ? eXsum_Raw/eWgt_Raw:0.0;
  1333. QVectorE_RawY = (eWgt_Raw > 0.0) ? eYsum_Raw/eWgt_Raw:0.0;
  1334. QEastRaw.Set(QVectorE_RawX,QVectorE_RawY);
  1335. qRawEast[iHarm] = QEastRaw;
  1336. if(qRawEast[iHarm].Mod()){
  1337. psiRawEastBBC[iHarm] = qRawEast[iHarm].Phi()/(iHarm+1.);
  1338. if(psiRawEastBBC[iHarm] < 0.0) {psiRawEastBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1339. }
  1340. }
  1341. for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
  1342. {
  1343. hBBCFullQxRaw[iHarm]->Fill(qRawFull[iHarm].X(),event->cent16());
  1344. hBBCFullQyRaw[iHarm]->Fill(qRawFull[iHarm].Y(),event->cent16());
  1345. hBBCEastQxRaw[iHarm]->Fill(qRawEast[iHarm].X(),event->cent16());
  1346. hBBCEastQyRaw[iHarm]->Fill(qRawEast[iHarm].Y(),event->cent16());
  1347. hBBCWestQxRaw[iHarm]->Fill(qRawWest[iHarm].X(),event->cent16());
  1348. hBBCWestQyRaw[iHarm]->Fill(qRawWest[iHarm].Y(),event->cent16());
  1349. hBBCFullPsiRaw[iHarm]->Fill(psiRawFullBBC[iHarm], event->cent16());
  1350. hBBCEastPsiRaw[iHarm]->Fill(psiRawEastBBC[iHarm], event->cent16());
  1351. hBBCWestPsiRaw[iHarm]->Fill(psiRawWestBBC[iHarm], event->cent16());
  1352. }
  1353. // Gain & Rec & Shift
  1354. Double_t psiGainEastBBC[fNharmonics+1], psiGainWestBBC[fNharmonics+1], psiGainFullBBC[fNharmonics+1];
  1355. Double_t psiRecEastBBC[fNharmonics+1], psiRecWestBBC[fNharmonics+1], psiRecFullBBC[fNharmonics+1];
  1356. Double_t adcNormEastInner = BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),1,6)/6.;
  1357. Double_t adcNormEastOuter = BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),7,16)/10.;
  1358. Double_t adcNormWestInner = BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),1,6)/6.;
  1359. Double_t adcNormWestOuter = BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),7,16)/10.;
  1360. float egain[16]={
  1361. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),1))/adcNormEastInner),
  1362. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),2))/adcNormEastInner),
  1363. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),3))/adcNormEastInner),
  1364. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),4))/adcNormEastInner),
  1365. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),5))/adcNormEastInner),
  1366. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),6))/adcNormEastInner),
  1367. (float)((0.5)*BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),7))/adcNormEastOuter),
  1368. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),8))/adcNormEastOuter),
  1369. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),9))/adcNormEastOuter),
  1370. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),10))/adcNormEastOuter),
  1371. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),11))/adcNormEastOuter),
  1372. (float)((0.5)*BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),12))/adcNormEastOuter),
  1373. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),13))/adcNormEastOuter),
  1374. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),14))/adcNormEastOuter),
  1375. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),15))/adcNormEastOuter),
  1376. (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),16))/adcNormEastOuter)
  1377. };
  1378. float wgain[16]={
  1379. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),1))/adcNormWestInner),
  1380. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),2))/adcNormWestInner),
  1381. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),3))/adcNormWestInner),
  1382. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),4))/adcNormWestInner),
  1383. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),5))/adcNormWestInner),
  1384. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),6))/adcNormWestInner),
  1385. (float)((0.5)*BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),7))/adcNormWestOuter),
  1386. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),8))/adcNormWestOuter),
  1387. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),9))/adcNormWestOuter),
  1388. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),10))/adcNormWestOuter),
  1389. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),11))/adcNormWestOuter),
  1390. (float)((0.5)*BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),12))/adcNormWestOuter),
  1391. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),13))/adcNormWestOuter),
  1392. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),14))/adcNormWestOuter),
  1393. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),15))/adcNormWestOuter),
  1394. (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),16))/adcNormWestOuter)
  1395. };
  1396. float nHitE_Gain[16], nHitW_Gain[16];
  1397. for(int i=0;i<16;i++){
  1398. nHitE_Gain[i] = event->bbcAdcEast(i)/egain[i];
  1399. nHitW_Gain[i] = event->bbcAdcWest(i)/wgain[i];
  1400. }
  1401. Double_t psiWest[fNharmonics+1];
  1402. Double_t psiEast[fNharmonics+1];
  1403. Double_t psiFull[fNharmonics+1];
  1404. TVector2 qGainEast[fNharmonics+1], qGainWest[fNharmonics+1], qGainFull[fNharmonics+1];
  1405. TVector2 qRecEast[fNharmonics+1], qRecWest[fNharmonics+1], qRecFull[fNharmonics+1];
  1406. for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
  1407. {
  1408. //BBC Full
  1409. //v1*y > 0 by conventions most reasonable to the STAR event plane
  1410. //therefore the west event plane has the right sign for the event plane
  1411. //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
  1412. //The east event plane points opposite that of the west and the west finds the "correct" EP for the STAR coordinates, use West-East
  1413. Double_t QVector_GainX = 0., QVector_GainY = 0.;
  1414. Double_t EventPlaneGainFull = 0.;
  1415. Float_t eXFsum_Gain = 0., wXFsum_Gain = 0., eYFsum_Gain = 0., wYFsum_Gain = 0., eWFgt_Gain=0., wWFgt_Gain = 0.;
  1416. TVector2 QFullGain;
  1417. for(int iTile = 0; iTile < 16; iTile++) {
  1418. eXFsum_Gain += cos((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
  1419. wXFsum_Gain += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
  1420. eYFsum_Gain += sin((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
  1421. wYFsum_Gain += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
  1422. eWFgt_Gain += nHitE_Gain[iTile];
  1423. wWFgt_Gain += nHitW_Gain[iTile];
  1424. }
  1425. QVector_GainX=(eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wXFsum_Gain/wWFgt_Gain + eXFsum_Gain/eWFgt_Gain:0.;
  1426. QVector_GainY=(eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wYFsum_Gain/wWFgt_Gain + eYFsum_Gain/eWFgt_Gain:0.;
  1427. QFullGain.Set(QVector_GainX,QVector_GainY);
  1428. if(QFullGain.Mod()){
  1429. EventPlaneGainFull = QFullGain.Phi();
  1430. if(EventPlaneGainFull < 0.0) {EventPlaneGainFull += 2.*TMath::Pi()/(iHarm+1.);}
  1431. }
  1432. psiFull[iHarm] = EventPlaneGainFull;
  1433. qGainFull[iHarm] = QFullGain;
  1434. if(qGainFull[iHarm].Mod()){
  1435. psiGainFullBBC[iHarm] = qGainFull[iHarm].Phi()/(iHarm+1.);
  1436. if(psiGainFullBBC[iHarm] < 0.0) {psiGainFullBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1437. }
  1438. //BBC West
  1439. //v1*y > 0 by conventions most reasonable to the STAR event plane
  1440. //therefore the west event plane has the right sign for the event plane
  1441. //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
  1442. Double_t QVectorW_GainX = 0., QVectorW_GainY = 0.;
  1443. Double_t EventPlaneGainWest = 0.;
  1444. Float_t wXsum_Gain = 0., wYsum_Gain = 0., wWgt_Gain = 0.;
  1445. TVector2 QWestGain;
  1446. for(int iTile = 0; iTile < 16; iTile++) {
  1447. wXsum_Gain += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
  1448. wYsum_Gain += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
  1449. wWgt_Gain += nHitW_Gain[iTile];
  1450. }
  1451. QVectorW_GainX = (wWgt_Gain > 0.0) ? wXsum_Gain/wWgt_Gain:0.0;
  1452. QVectorW_GainY = (wWgt_Gain > 0.0) ? wYsum_Gain/wWgt_Gain:0.0;
  1453. QWestGain.Set(QVectorW_GainX,QVectorW_GainY);
  1454. if(QWestGain.Mod()){
  1455. EventPlaneGainWest = QWestGain.Phi()/(iHarm+1.);
  1456. if(EventPlaneGainWest < 0.0) {EventPlaneGainWest += 2.*TMath::Pi()/(iHarm+1.);}
  1457. }
  1458. psiWest[iHarm] = EventPlaneGainWest;
  1459. qGainWest[iHarm] = QWestGain;
  1460. if(qGainWest[iHarm].Mod()){
  1461. psiGainWestBBC[iHarm] = qGainWest[iHarm].Phi()/(iHarm+1.);
  1462. if(psiGainWestBBC[iHarm] < 0.0) {psiGainWestBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1463. }
  1464. //BBC East
  1465. //in BBC_GetPhi the fact that the tile numbering scheme is reversed about the y axis for the east BBC is accounted for
  1466. //The east event plane points opposite that of the west since the spectators are on different sides -- see west comments
  1467. Double_t QVectorE_GainX = 0., QVectorE_GainY = 0.;
  1468. Double_t EventPlaneGainEast = 0.;
  1469. Float_t eXsum_Gain = 0., eYsum_Gain = 0., eWgt_Gain = 0.;
  1470. TVector2 QEastGain;
  1471. for(int iTile = 0; iTile < 16; iTile++) {
  1472. eXsum_Gain += cos((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
  1473. eYsum_Gain += sin((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
  1474. eWgt_Gain += nHitE_Gain[iTile];
  1475. }
  1476. QVectorE_GainX = (eWgt_Gain > 0.0) ? eXsum_Gain/eWgt_Gain:0.0;
  1477. QVectorE_GainY = (eWgt_Gain > 0.0) ? eYsum_Gain/eWgt_Gain:0.0;
  1478. QEastGain.Set(QVectorE_GainX,QVectorE_GainY);
  1479. if(QEastGain.Mod()){
  1480. EventPlaneGainEast = QEastGain.Phi()/(iHarm+1.);
  1481. if(EventPlaneGainEast < 0.0) {EventPlaneGainEast += 2.*TMath::Pi()/(iHarm+1.);}
  1482. }
  1483. psiEast[iHarm] = EventPlaneGainEast;
  1484. qGainEast[iHarm] = QEastGain;
  1485. if(qGainEast[iHarm].Mod()){
  1486. psiGainEastBBC[iHarm] = qGainEast[iHarm].Phi()/(iHarm+1.);
  1487. if(psiGainEastBBC[iHarm] < 0.0) {psiGainEastBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1488. }
  1489. }
  1490. // Apply recentering
  1491. for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
  1492. {
  1493. TVector2 recFullTMP((Double_t)p_BBC_Qx_Full_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
  1494. (Double_t)p_BBC_Qy_Full_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1495. TVector2 recEastTMP((Double_t)p_BBC_Qx_East_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
  1496. (Double_t)p_BBC_Qy_East_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1497. TVector2 recWestTMP((Double_t)p_BBC_Qx_West_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
  1498. (Double_t)p_BBC_Qy_West_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1499. qRecFull[iHarm] = qGainFull[iHarm] - recFullTMP;
  1500. qRecEast[iHarm] = qGainEast[iHarm] - recEastTMP;
  1501. qRecWest[iHarm] = qGainWest[iHarm] - recWestTMP;
  1502. if(qRecFull[iHarm].Mod()){
  1503. psiRecFullBBC[iHarm] = qRecFull[iHarm].Phi()/(iHarm+1.);
  1504. if(psiRecFullBBC[iHarm] < 0.0) {psiRecFullBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1505. }
  1506. if(qRecEast[iHarm].Mod()){
  1507. psiRecEastBBC[iHarm] = qRecEast[iHarm].Phi()/(iHarm+1.);
  1508. if(psiRecEastBBC[iHarm] < 0.0) {psiRecEastBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1509. }
  1510. if(qRecWest[iHarm].Mod()){
  1511. psiRecWestBBC[iHarm] = qRecWest[iHarm].Phi()/(iHarm+1.);
  1512. if(psiRecWestBBC[iHarm] < 0.0) {psiRecWestBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
  1513. }
  1514. }
  1515. for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
  1516. {
  1517. // Full BBC
  1518. if (qRecFull[iHarm].Mod())
  1519. {
  1520. psiEP_rec = qRecFull[iHarm].Phi() / (iHarm+1.);
  1521. if (psiEP_rec < 0.) {psiEP_rec += 2*TMath::Pi()/(iHarm+1.);}
  1522. }
  1523. dPsi = 0.;
  1524. for (int iK = 0; iK < NShiftOrderMax*3; iK++)
  1525. {
  1526. shiftCos = p_BBC_Cos_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1527. shiftSin = p_BBC_Sin_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1528. dPsi += (1. / (iHarm+1.)) * (2. / (iK + 1.)) *
  1529. (-1. * shiftSin * cos((iHarm+1.) * (iK + 1.) * psiEP_rec) +
  1530. 1. * shiftCos * sin((iHarm+1.) * (iK + 1.) * psiEP_rec));
  1531. }
  1532. psiEP_BBC_Full[iHarm] = psiEP_rec + dPsi;
  1533. // East BBC
  1534. if (qRecEast[iHarm].Mod())
  1535. {
  1536. psiEP_rec = qRecEast[iHarm].Phi() / (iHarm+1.);
  1537. if (psiEP_rec < 0.) {psiEP_rec += 2*TMath::Pi()/(iHarm+1.);}
  1538. }
  1539. dPsi = 0.;
  1540. for (int iK = 0; iK < NShiftOrderMax*3; iK++)
  1541. {
  1542. shiftCos = p_BBC_Cos_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1543. shiftSin = p_BBC_Sin_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1544. dPsi += (1. / (iHarm+1.)) * (2. / (iK + 1.)) *
  1545. (-1. * shiftSin * cos((iHarm+1.) * (iK + 1.) * psiEP_rec) +
  1546. 1. * shiftCos * sin((iHarm+1.) * (iK + 1.) * psiEP_rec));
  1547. }
  1548. psiEP_BBC_East[iHarm] = psiEP_rec + dPsi;
  1549. // West BBC
  1550. if (qRecWest[iHarm].Mod())
  1551. {
  1552. psiEP_rec = qRecWest[iHarm].Phi() / (iHarm+1.);
  1553. if (psiEP_rec < 0.) {psiEP_rec += 2*TMath::Pi()/(iHarm+1.);}
  1554. }
  1555. dPsi = 0.;
  1556. for (int iK = 0; iK < NShiftOrderMax*3; iK++)
  1557. {
  1558. shiftCos = p_BBC_Cos_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1559. shiftSin = p_BBC_Sin_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1560. dPsi += (1. / (iHarm+1.)) * (2. / (iK + 1.)) *
  1561. (-1. * shiftSin * cos((iHarm+1.) * (iK + 1.) * psiEP_rec) +
  1562. 1. * shiftCos * sin((iHarm+1.) * (iK + 1.) * psiEP_rec));
  1563. }
  1564. psiEP_BBC_West[iHarm] = psiEP_rec + dPsi;
  1565. if(psiEP_BBC_Full[iHarm] < 0.0) psiEP_BBC_Full[iHarm] += 2.*TMath::Pi()/(iHarm+1.);
  1566. if(psiEP_BBC_Full[iHarm] > 2.*TMath::Pi()/(iHarm+1.)) psiEP_BBC_Full[iHarm] -= 2.*TMath::Pi()/(iHarm+1.);
  1567. if(psiEP_BBC_East[iHarm] < 0.0) psiEP_BBC_East[iHarm] += 2.*TMath::Pi()/(iHarm+1.);
  1568. if(psiEP_BBC_East[iHarm] > 2.*TMath::Pi()/(iHarm+1.)) psiEP_BBC_East[iHarm] -= 2.*TMath::Pi()/(iHarm+1.);
  1569. if(psiEP_BBC_West[iHarm] < 0.0) psiEP_BBC_West[iHarm] += 2.*TMath::Pi()/(iHarm+1.);
  1570. if(psiEP_BBC_West[iHarm] > 2.*TMath::Pi()/(iHarm+1.)) psiEP_BBC_West[iHarm] -= 2.*TMath::Pi()/(iHarm+1.);
  1571. }
  1572. for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
  1573. {
  1574. hBBCFullQxGain[iHarm]->Fill(qGainFull[iHarm].X(),event->cent16());
  1575. hBBCFullQyGain[iHarm]->Fill(qGainFull[iHarm].Y(),event->cent16());
  1576. hBBCEastQxGain[iHarm]->Fill(qGainEast[iHarm].X(),event->cent16());
  1577. hBBCEastQyGain[iHarm]->Fill(qGainEast[iHarm].Y(),event->cent16());
  1578. hBBCWestQxGain[iHarm]->Fill(qGainWest[iHarm].X(),event->cent16());
  1579. hBBCWestQyGain[iHarm]->Fill(qGainWest[iHarm].Y(),event->cent16());
  1580. hBBCFullPsiGain[iHarm]->Fill(psiGainFullBBC[iHarm], event->cent16());
  1581. hBBCEastPsiGain[iHarm]->Fill(psiGainEastBBC[iHarm], event->cent16());
  1582. hBBCWestPsiGain[iHarm]->Fill(psiGainWestBBC[iHarm], event->cent16());
  1583. hBBCFullQxRec[iHarm]->Fill(qRecFull[iHarm].X(),event->cent16());
  1584. hBBCFullQyRec[iHarm]->Fill(qRecFull[iHarm].Y(),event->cent16());
  1585. hBBCEastQxRec[iHarm]->Fill(qRecEast[iHarm].X(),event->cent16());
  1586. hBBCEastQyRec[iHarm]->Fill(qRecEast[iHarm].Y(),event->cent16());
  1587. hBBCWestQxRec[iHarm]->Fill(qRecWest[iHarm].X(),event->cent16());
  1588. hBBCWestQyRec[iHarm]->Fill(qRecWest[iHarm].Y(),event->cent16());
  1589. hBBCFullPsiRec[iHarm]->Fill(psiRecFullBBC[iHarm], event->cent16());
  1590. hBBCEastPsiRec[iHarm]->Fill(psiRecEastBBC[iHarm], event->cent16());
  1591. hBBCWestPsiRec[iHarm]->Fill(psiRecWestBBC[iHarm], event->cent16());
  1592. hBBCFullPsiShift[iHarm]->Fill(psiEP_BBC_Full[iHarm], event->cent16());
  1593. hBBCEastPsiShift[iHarm]->Fill(psiEP_BBC_East[iHarm], event->cent16());
  1594. hBBCWestPsiShift[iHarm]->Fill(psiEP_BBC_West[iHarm], event->cent16());
  1595. // East
  1596. pQx2BBCvsRun[0]->Fill(event->runId(),qRecEast[1].X());
  1597. pQy2BBCvsRun[0]->Fill(event->runId(),qRecEast[1].Y());
  1598. pCosBBCvsRun[0][0]->Fill(event->runId(),TMath::Cos(1*psiEP_BBC_East[1]));
  1599. pCosBBCvsRun[1][0]->Fill(event->runId(),TMath::Cos(2*psiEP_BBC_East[1]));
  1600. pCosBBCvsRun[2][0]->Fill(event->runId(),TMath::Cos(3*psiEP_BBC_East[1]));
  1601. pCosBBCvsRun[3][0]->Fill(event->runId(),TMath::Cos(4*psiEP_BBC_East[1]));
  1602. pSinBBCvsRun[0][0]->Fill(event->runId(),TMath::Sin(1*psiEP_BBC_East[1]));
  1603. pSinBBCvsRun[1][0]->Fill(event->runId(),TMath::Sin(2*psiEP_BBC_East[1]));
  1604. pSinBBCvsRun[2][0]->Fill(event->runId(),TMath::Sin(3*psiEP_BBC_East[1]));
  1605. pSinBBCvsRun[3][0]->Fill(event->runId(),TMath::Sin(4*psiEP_BBC_East[1]));
  1606. // West
  1607. pQx2BBCvsRun[1]->Fill(event->runId(),qRecWest[1].X());
  1608. pQy2BBCvsRun[1]->Fill(event->runId(),qRecWest[1].Y());
  1609. pCosBBCvsRun[0][1]->Fill(event->runId(),TMath::Cos(1*psiEP_BBC_West[1]));
  1610. pCosBBCvsRun[1][1]->Fill(event->runId(),TMath::Cos(2*psiEP_BBC_West[1]));
  1611. pCosBBCvsRun[2][1]->Fill(event->runId(),TMath::Cos(3*psiEP_BBC_West[1]));
  1612. pCosBBCvsRun[3][1]->Fill(event->runId(),TMath::Cos(4*psiEP_BBC_West[1]));
  1613. pSinBBCvsRun[0][1]->Fill(event->runId(),TMath::Sin(1*psiEP_BBC_West[1]));
  1614. pSinBBCvsRun[1][1]->Fill(event->runId(),TMath::Sin(2*psiEP_BBC_West[1]));
  1615. pSinBBCvsRun[2][1]->Fill(event->runId(),TMath::Sin(3*psiEP_BBC_West[1]));
  1616. pSinBBCvsRun[3][1]->Fill(event->runId(),TMath::Sin(4*psiEP_BBC_West[1]));
  1617. }
  1618. // Track analysis
  1619. Int_t nTracks = dst->numberOfTracks();
  1620. // Track loop
  1621. for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
  1622. {
  1623. // Retrieve i-th femto track
  1624. StFemtoTrack *femtoTrack = dst->track(iTrk);
  1625. if (!femtoTrack)
  1626. continue;
  1627. hNhitsFit[0]->Fill(femtoTrack->nHitsFit());
  1628. hNhitsPoss[0]->Fill(femtoTrack->nHitsPoss());
  1629. if (femtoTrack->nHitsPoss() != 0)
  1630. {
  1631. hNhitsRatio[0]->Fill((Double_t)(femtoTrack->nHitsFit() / femtoTrack->nHitsPoss()));
  1632. }
  1633. hPt[0]->Fill(femtoTrack->pt());
  1634. hPtot[0]->Fill(femtoTrack->ptot());
  1635. hEta[0]->Fill(femtoTrack->eta());
  1636. hPhi[0]->Fill(femtoTrack->phi());
  1637. hDCA[0]->Fill(femtoTrack->gDCA(pVtx).Mag());
  1638. pPtvsRun[0]->Fill(event->runId(),femtoTrack->pt());
  1639. pEtavsRun[0]->Fill(event->runId(),femtoTrack->eta());
  1640. pPhivsRun[0]->Fill(event->runId(),femtoTrack->phi());
  1641. pNhitsFitvsRun[0]->Fill(event->runId(),femtoTrack->nHitsFit());
  1642. pNhitsPossvsRun[0]->Fill(event->runId(),femtoTrack->nHitsPoss());
  1643. pDCAvsRun[0]->Fill(event->runId(),femtoTrack->gDCA(pVtx).Mag());
  1644. pdEdxvsRun[0]->Fill(event->runId(),femtoTrack->dEdx());
  1645. // Must be a primary track
  1646. if (!femtoTrack->isPrimary())
  1647. continue;
  1648. //Track selection
  1649. if (isGoodTrack(femtoTrack, energy, pVtx))
  1650. {
  1651. fUsedTracks++;
  1652. //Fill recentered Q-vectors for all TPC (FULL)
  1653. weight = GetWeight(femtoTrack);
  1654. QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1655. QvTMP.SetY(p_Qx_FULL_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1656. fQv2FullEP += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
  1657. QvTMP.Set(0., 0.);
  1658. QvTMP.SetX(p_Qx_FULL_EP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1659. QvTMP.SetY(p_Qy_FULL_EP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1660. fQv3FullEP += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
  1661. QvTMP.Set(0., 0.);
  1662. QvTMP.SetX(p_Qx_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1663. QvTMP.SetY(p_Qy_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1664. fQv2FullSP += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
  1665. QvTMP.Set(0., 0.);
  1666. QvTMP.SetX(p_Qx_FULL_SP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1667. QvTMP.SetY(p_Qy_FULL_SP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
  1668. fQv3FullSP += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
  1669. QvTMP.Set(0.,0.);
  1670. fQcountFull++;
  1671. if (femtoTrack->eta() >= 0.)
  1672. {
  1673. fQcountFullWest++;
  1674. }
  1675. else
  1676. {
  1677. fQcountFullEast++;
  1678. }
  1679. for (int iGap=0; iGap<NEtaGaps; iGap++)
  1680. {
  1681. //fill Qvectors from EAST arm
  1682. if (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < 0.)//cutEtaGap[iGap])
  1683. {
  1684. QvTMP.SetX(p_Qx_EAST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1685. QvTMP.SetY(p_Qy_EAST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1686. fQv2EastEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
  1687. QvTMP.Set(0., 0.);
  1688. QvTMP.SetX(p_Qx_EAST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1689. QvTMP.SetY(p_Qy_EAST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1690. fQv3EastEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
  1691. QvTMP.Set(0., 0.);
  1692. QvTMP.SetX(p_Qx_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1693. QvTMP.SetY(p_Qy_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1694. fQv2EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
  1695. QvTMP.Set(0., 0.);
  1696. QvTMP.SetX(p_Qx_EAST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1697. QvTMP.SetY(p_Qy_EAST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1698. fQv3EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
  1699. QvTMP.Set(0.,0.);
  1700. fQcountEast[iGap]++;
  1701. }
  1702. //fill Qvectors from WEST arm
  1703. if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 0.)//cutEtaGap[iGap])
  1704. {
  1705. QvTMP.SetX(p_Qx_WEST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1706. QvTMP.SetY(p_Qy_WEST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1707. fQv2WestEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
  1708. QvTMP.Set(0., 0.);
  1709. QvTMP.SetX(p_Qx_WEST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1710. QvTMP.SetY(p_Qy_WEST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1711. fQv3WestEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
  1712. QvTMP.Set(0., 0.);
  1713. QvTMP.SetX(p_Qx_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1714. QvTMP.SetY(p_Qy_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1715. fQv2WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
  1716. QvTMP.Set(0., 0.);
  1717. QvTMP.SetX(p_Qx_WEST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1718. QvTMP.SetY(p_Qy_WEST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
  1719. fQv3WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
  1720. QvTMP.Set(0.,0.);
  1721. fQcountWest[iGap]++;
  1722. }
  1723. }
  1724. }
  1725. //Track selection
  1726. if (!isGoodTrackFlow(femtoTrack, energy, pVtx))
  1727. continue;
  1728. hNhitsFit[1]->Fill(femtoTrack->nHitsFit());
  1729. hNhitsPoss[1]->Fill(femtoTrack->nHitsPoss());
  1730. if (femtoTrack->nHitsPoss() != 0)
  1731. {
  1732. hNhitsRatio[1]->Fill((Double_t)(femtoTrack->nHitsFit() / femtoTrack->nHitsPoss()));
  1733. }
  1734. hPt[1]->Fill(femtoTrack->pt());
  1735. hPtot[1]->Fill(femtoTrack->ptot());
  1736. hPhi[1]->Fill(femtoTrack->phi());
  1737. hEta[1]->Fill(femtoTrack->eta());
  1738. hDCA[1]->Fill(femtoTrack->gDCA(pVtx).Mag());
  1739. pPtvsRun[1]->Fill(event->runId(),femtoTrack->pt());
  1740. pEtavsRun[1]->Fill(event->runId(),femtoTrack->eta());
  1741. pPhivsRun[1]->Fill(event->runId(),femtoTrack->phi());
  1742. pNhitsFitvsRun[1]->Fill(event->runId(),femtoTrack->nHitsFit());
  1743. pNhitsPossvsRun[1]->Fill(event->runId(),femtoTrack->nHitsPoss());
  1744. pDCAvsRun[1]->Fill(event->runId(),femtoTrack->gDCA(pVtx).Mag());
  1745. pdEdxvsRun[1]->Fill(event->runId(),femtoTrack->dEdx());
  1746. Int_t i_part = -1;
  1747. //TPC-only
  1748. if (!femtoTrack->isTofTrack())
  1749. {
  1750. //pion id
  1751. if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.6 &&
  1752. TMath::Abs(femtoTrack->nSigmaPion()) < 2)
  1753. {
  1754. i_part = 0;
  1755. }
  1756. // kaon id
  1757. if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.5 &&
  1758. TMath::Abs(femtoTrack->nSigmaKaon()) < 2)
  1759. {
  1760. i_part = 1;
  1761. }
  1762. // proton id
  1763. if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 0.9 &&
  1764. TMath::Abs(femtoTrack->nSigmaProton()) < 2)
  1765. {
  1766. i_part = 2;
  1767. }
  1768. //pion id
  1769. if (femtoTrack->ptot() >= 0.6 && femtoTrack->ptot() < 0.7 &&
  1770. TMath::Abs(femtoTrack->nSigmaPion()) < 2 &&
  1771. TMath::Abs(femtoTrack->nSigmaKaon()) > 2)
  1772. {
  1773. i_part = 0;
  1774. }
  1775. // kaon id
  1776. if (femtoTrack->ptot() >= 0.5 && femtoTrack->ptot() < 0.7 &&
  1777. TMath::Abs(femtoTrack->nSigmaKaon()) < 2 &&
  1778. TMath::Abs(femtoTrack->nSigmaPion()) > 3)
  1779. {
  1780. i_part = 1;
  1781. }
  1782. // proton id
  1783. if (femtoTrack->ptot() >= 0.9 && femtoTrack->ptot() < 1.2 &&
  1784. TMath::Abs(femtoTrack->nSigmaProton()) < 2 &&
  1785. TMath::Abs(femtoTrack->nSigmaPion()) > 3)
  1786. {
  1787. i_part = 2;
  1788. }
  1789. }
  1790. //TPC+TOF
  1791. if (isGoodPID(femtoTrack))
  1792. {
  1793. // pion id
  1794. if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
  1795. TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
  1796. femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)
  1797. {
  1798. i_part = 0;
  1799. }
  1800. // kaon id
  1801. if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
  1802. TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
  1803. femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)
  1804. {
  1805. i_part = 1;
  1806. }
  1807. // proton id
  1808. if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
  1809. TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
  1810. femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)
  1811. {
  1812. i_part = 2;
  1813. }
  1814. }
  1815. for (int iGap=0; iGap<NEtaGaps; iGap++)
  1816. {
  1817. //EAST OR WEST
  1818. if ( (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < -1.*cutEtaGap.at(iGap)) ||
  1819. (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 1.*cutEtaGap.at(iGap)) )
  1820. {
  1821. //Fill recentered Q-vectors for all TPC (FULL)
  1822. weight = GetWeight(femtoTrack);
  1823. h_pt_spectra[3][(Int_t)event->cent16()]->Fill((Double_t)event->runId(),femtoTrack->pt(),weight);
  1824. if (i_part != -1)
  1825. {
  1826. h_pt_spectra[i_part][(Int_t)event->cent16()]->Fill((Double_t)event->runId(),femtoTrack->pt(),weight);
  1827. }
  1828. }
  1829. }
  1830. // Check if track has TOF signal
  1831. if (femtoTrack->isTofTrack())
  1832. {
  1833. if (femtoTrack->gDCA(pVtx).Mag() >= 1.)
  1834. continue;
  1835. if (femtoTrack->pt() > 0.5 && femtoTrack->pt() < 0.7)
  1836. {
  1837. hnSigM2_lowpt[0]->Fill(femtoTrack->nSigmaPion(), femtoTrack->massSqr());
  1838. }
  1839. if (femtoTrack->pt() > 2.2 && femtoTrack->pt() < 2.4)
  1840. {
  1841. hnSigM2_highpt[0]->Fill(femtoTrack->nSigmaPion(), femtoTrack->massSqr());
  1842. }
  1843. hdEdxQp[0][0]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
  1844. hM2Qp[0][0]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
  1845. hNsigQp[0][0]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
  1846. hM2pt[0][0]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
  1847. hNsigpt[0][0]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
  1848. // pion
  1849. if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
  1850. TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
  1851. femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)
  1852. {
  1853. hdEdxQp[0][1]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
  1854. hM2Qp[0][1]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
  1855. hNsigQp[0][1]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
  1856. hM2pt[0][1]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
  1857. hNsigpt[0][1]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
  1858. }
  1859. // kaon
  1860. if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
  1861. TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
  1862. femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)
  1863. {
  1864. hdEdxQp[0][2]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
  1865. hM2Qp[0][2]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
  1866. hNsigQp[0][2]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
  1867. hM2pt[0][2]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
  1868. hNsigpt[0][2]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
  1869. }
  1870. // proton
  1871. if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
  1872. TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
  1873. femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)
  1874. {
  1875. hdEdxQp[0][3]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
  1876. hM2Qp[0][3]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
  1877. hNsigQp[0][3]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
  1878. hM2pt[0][3]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
  1879. hNsigpt[0][3]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
  1880. }
  1881. } //if( isTofTrack() )
  1882. } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
  1883. //Getting PsiEP with recentering + shift correction
  1884. Psi_ReCenter=0.;
  1885. Double_t mCos, mSin;
  1886. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  1887. {
  1888. Psi2FullEP = 0.;
  1889. Psi3FullEP = 0.;
  1890. //FULL 2-nd order
  1891. Psi_ReCenter = TMath::ATan2(fQv2FullEP.Y(), fQv2FullEP.X()) / 2.;
  1892. dPsi = 0.;
  1893. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1894. {
  1895. mCos = p_Cos2_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1896. mSin = p_Cos2_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1897. dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
  1898. }
  1899. Psi2FullEP = GetEPAngle(Psi_ReCenter + dPsi, 2.);
  1900. //FULL 3-nd order
  1901. Psi_ReCenter = TMath::ATan2(fQv3FullEP.Y(), fQv3FullEP.X()) / 3.;
  1902. dPsi = 0.;
  1903. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1904. {
  1905. mCos = p_Cos3_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1906. mSin = p_Cos3_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
  1907. dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
  1908. }
  1909. Psi3FullEP = GetEPAngle(Psi_ReCenter + dPsi, 3.);
  1910. }
  1911. for (int iGap=0; iGap<NEtaGaps; iGap++)
  1912. {
  1913. if (fQcountEast[iGap] > cutNcountQvGap &&
  1914. fQcountWest[iGap] > cutNcountQvGap)
  1915. {
  1916. Psi2EastEP[iGap]=0.;
  1917. Psi3EastEP[iGap]=0.;
  1918. Psi2WestEP[iGap]=0.;
  1919. Psi3WestEP[iGap]=0.;
  1920. Psi2EastSP[iGap]=0.;
  1921. Psi3EastSP[iGap]=0.;
  1922. Psi2WestSP[iGap]=0.;
  1923. Psi3WestSP[iGap]=0.;
  1924. //EAST 2-nd order
  1925. Psi_ReCenter = TMath::ATan2(fQv2EastEP[iGap].Y(),fQv2EastEP[iGap].X())/2.;
  1926. dPsi = 0.;
  1927. for (int iK=0; iK<NShiftOrderMax; iK++)
  1928. {
  1929. mCos = p_Cos2_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1930. mSin = p_Cos2_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1931. dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
  1932. TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
  1933. mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
  1934. }
  1935. Psi2EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
  1936. Psi_ReCenter = TMath::ATan2(fQv2EastSP[iGap].Y(), fQv2EastSP[iGap].X()) / 2.;
  1937. dPsi = 0.;
  1938. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1939. {
  1940. mCos = p_Cos2_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1941. mSin = p_Cos2_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1942. dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
  1943. }
  1944. Psi2EastSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
  1945. fQv2EastSP_Shift[iGap].Set(fQv2EastSP[iGap].Mod()*TMath::Cos(2.*Psi2EastSP[iGap]),fQv2EastSP[iGap].Mod()*TMath::Sin(2.*Psi2EastSP[iGap]));
  1946. //WEST 2-nd order
  1947. Psi_ReCenter = TMath::ATan2(fQv2WestEP[iGap].Y(),fQv2WestEP[iGap].X())/2.;
  1948. dPsi = 0.;
  1949. for (int iK=0; iK<NShiftOrderMax; iK++)
  1950. {
  1951. mCos = p_Cos2_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1952. mSin = p_Cos2_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1953. dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
  1954. TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
  1955. mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
  1956. }
  1957. Psi2WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
  1958. Psi_ReCenter = TMath::ATan2(fQv2WestSP[iGap].Y(), fQv2WestSP[iGap].X()) / 2.;
  1959. dPsi = 0.;
  1960. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1961. {
  1962. mCos = p_Cos2_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1963. mSin = p_Cos2_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1964. dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
  1965. }
  1966. Psi2WestSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
  1967. fQv2WestSP_Shift[iGap].Set(fQv2WestSP[iGap].Mod()*TMath::Cos(2.*Psi2WestSP[iGap]),fQv2WestSP[iGap].Mod()*TMath::Sin(2.*Psi2WestSP[iGap]));
  1968. //EAST 3-nd order
  1969. Psi_ReCenter = TMath::ATan2(fQv3EastEP[iGap].Y(),fQv3EastEP[iGap].X())/3.;
  1970. dPsi = 0.;
  1971. for (int iK=0; iK<NShiftOrderMax; iK++)
  1972. {
  1973. mCos = p_Cos3_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1974. mSin = p_Cos3_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1975. dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
  1976. TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
  1977. mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
  1978. }
  1979. Psi3EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
  1980. Psi_ReCenter = TMath::ATan2(fQv3EastSP[iGap].Y(), fQv3EastSP[iGap].X()) / 3.;
  1981. dPsi = 0.;
  1982. for (int iK = 0; iK < NShiftOrderMax; iK++)
  1983. {
  1984. mCos = p_Cos3_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1985. mSin = p_Cos3_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1986. dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
  1987. }
  1988. Psi3EastSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
  1989. fQv3EastSP_Shift[iGap].Set(fQv3EastSP[iGap].Mod()*TMath::Cos(3.*Psi3EastSP[iGap]),fQv3EastSP[iGap].Mod()*TMath::Sin(3.*Psi3EastSP[iGap]));
  1990. //WEST 3-nd order
  1991. Psi_ReCenter = TMath::ATan2(fQv3WestEP[iGap].Y(),fQv3WestEP[iGap].X())/3.;
  1992. dPsi = 0.;
  1993. for (int iK=0; iK<NShiftOrderMax; iK++)
  1994. {
  1995. mCos = p_Cos3_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1996. mSin = p_Cos3_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  1997. dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
  1998. TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
  1999. mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
  2000. }
  2001. Psi3WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
  2002. Psi_ReCenter = TMath::ATan2(fQv3WestSP[iGap].Y(), fQv3WestSP[iGap].X()) / 3.;
  2003. dPsi = 0.;
  2004. for (int iK = 0; iK < NShiftOrderMax; iK++)
  2005. {
  2006. mCos = p_Cos3_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  2007. mSin = p_Cos3_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
  2008. dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
  2009. }
  2010. Psi3WestSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
  2011. fQv3WestSP_Shift[iGap].Set(fQv3WestSP[iGap].Mod()*TMath::Cos(3.*Psi3WestSP[iGap]),fQv3WestSP[iGap].Mod()*TMath::Sin(3.*Psi3WestSP[iGap]));
  2012. }
  2013. }
  2014. // 2nd harmonics TPC
  2015. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  2016. {
  2017. hTPCFullQxRaw[0]->Fill(fQv2FullEPRaw.X(),event->cent16());
  2018. hTPCFullQyRaw[0]->Fill(fQv2FullEPRaw.X(),event->cent16());
  2019. hTPCFullPsiRaw[0]->Fill(GetEPAngle(fQv2FullEPRaw,2), event->cent16());
  2020. }
  2021. if (fQcountEast[0] > cutNcountQvGap &&
  2022. fQcountWest[0] > cutNcountQvGap)
  2023. {
  2024. hTPCEastQxRaw[0]->Fill(fQv2EastEPRaw[0].X(),event->cent16());
  2025. hTPCEastQyRaw[0]->Fill(fQv2EastEPRaw[0].Y(),event->cent16());
  2026. hTPCWestQxRaw[0]->Fill(fQv2WestEPRaw[0].X(),event->cent16());
  2027. hTPCWestQyRaw[0]->Fill(fQv2WestEPRaw[0].Y(),event->cent16());
  2028. hTPCEastPsiRaw[0]->Fill(GetEPAngle(fQv2EastEPRaw[0],2), event->cent16());
  2029. hTPCWestPsiRaw[0]->Fill(GetEPAngle(fQv2WestEPRaw[0],2), event->cent16());
  2030. }
  2031. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  2032. {
  2033. hTPCFullQxRec[0]->Fill(fQv2FullEP.X(),event->cent16());
  2034. hTPCFullQyRec[0]->Fill(fQv2FullEP.X(),event->cent16());
  2035. hTPCFullPsiRec[0]->Fill(GetEPAngle(fQv2FullEP,2), event->cent16());
  2036. }
  2037. if (fQcountEast[0] > cutNcountQvGap &&
  2038. fQcountWest[0] > cutNcountQvGap)
  2039. {
  2040. hTPCEastQxRec[0]->Fill(fQv2EastEP[0].X(),event->cent16());
  2041. hTPCEastQyRec[0]->Fill(fQv2EastEP[0].Y(),event->cent16());
  2042. hTPCWestQxRec[0]->Fill(fQv2WestEP[0].X(),event->cent16());
  2043. hTPCWestQyRec[0]->Fill(fQv2WestEP[0].Y(),event->cent16());
  2044. hTPCEastPsiRec[0]->Fill(GetEPAngle(fQv2EastEP[0],2), event->cent16());
  2045. hTPCWestPsiRec[0]->Fill(GetEPAngle(fQv2WestEP[0],2), event->cent16());
  2046. }
  2047. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  2048. {
  2049. hTPCFullPsiShift[0]->Fill(GetEPAngle(Psi2FullEP,2), event->cent16());
  2050. // East
  2051. pQx2TPCvsRun[0]->Fill(event->runId(),fQv2EastEP[0].X());
  2052. pQy2TPCvsRun[0]->Fill(event->runId(),fQv2EastEP[0].Y());
  2053. pQx3TPCvsRun[0]->Fill(event->runId(),fQv3EastEP[0].X());
  2054. pQy3TPCvsRun[0]->Fill(event->runId(),fQv3EastEP[0].Y());
  2055. pCos2TPCvsRun[0][0]->Fill(event->runId(),TMath::Cos(1*(Psi2EastEP[0])));
  2056. pCos2TPCvsRun[1][0]->Fill(event->runId(),TMath::Cos(2*(Psi2EastEP[0])));
  2057. pCos2TPCvsRun[2][0]->Fill(event->runId(),TMath::Cos(3*(Psi2EastEP[0])));
  2058. pCos2TPCvsRun[3][0]->Fill(event->runId(),TMath::Cos(4*(Psi2EastEP[0])));
  2059. pCos3TPCvsRun[0][0]->Fill(event->runId(),TMath::Cos(1*(Psi3EastEP[0])));
  2060. pCos3TPCvsRun[1][0]->Fill(event->runId(),TMath::Cos(2*(Psi3EastEP[0])));
  2061. pCos3TPCvsRun[2][0]->Fill(event->runId(),TMath::Cos(3*(Psi3EastEP[0])));
  2062. pCos3TPCvsRun[3][0]->Fill(event->runId(),TMath::Cos(4*(Psi3EastEP[0])));
  2063. pSin2TPCvsRun[0][0]->Fill(event->runId(),TMath::Sin(1*(Psi2EastEP[0])));
  2064. pSin2TPCvsRun[1][0]->Fill(event->runId(),TMath::Sin(2*(Psi2EastEP[0])));
  2065. pSin2TPCvsRun[2][0]->Fill(event->runId(),TMath::Sin(3*(Psi2EastEP[0])));
  2066. pSin2TPCvsRun[3][0]->Fill(event->runId(),TMath::Sin(4*(Psi2EastEP[0])));
  2067. pSin3TPCvsRun[0][0]->Fill(event->runId(),TMath::Sin(1*(Psi3EastEP[0])));
  2068. pSin3TPCvsRun[1][0]->Fill(event->runId(),TMath::Sin(2*(Psi3EastEP[0])));
  2069. pSin3TPCvsRun[2][0]->Fill(event->runId(),TMath::Sin(3*(Psi3EastEP[0])));
  2070. pSin3TPCvsRun[3][0]->Fill(event->runId(),TMath::Sin(4*(Psi3EastEP[0])));
  2071. // West
  2072. pQx2TPCvsRun[1]->Fill(event->runId(),fQv2WestEP[0].X());
  2073. pQy2TPCvsRun[1]->Fill(event->runId(),fQv2WestEP[0].Y());
  2074. pQx3TPCvsRun[1]->Fill(event->runId(),fQv3WestEP[0].X());
  2075. pQy3TPCvsRun[1]->Fill(event->runId(),fQv3WestEP[0].Y());
  2076. pCos2TPCvsRun[0][1]->Fill(event->runId(),TMath::Cos(1*(Psi2WestEP[0])));
  2077. pCos2TPCvsRun[1][1]->Fill(event->runId(),TMath::Cos(2*(Psi2WestEP[0])));
  2078. pCos2TPCvsRun[2][1]->Fill(event->runId(),TMath::Cos(3*(Psi2WestEP[0])));
  2079. pCos2TPCvsRun[3][1]->Fill(event->runId(),TMath::Cos(4*(Psi2WestEP[0])));
  2080. pCos3TPCvsRun[0][1]->Fill(event->runId(),TMath::Cos(1*(Psi3WestEP[0])));
  2081. pCos3TPCvsRun[1][1]->Fill(event->runId(),TMath::Cos(2*(Psi3WestEP[0])));
  2082. pCos3TPCvsRun[2][1]->Fill(event->runId(),TMath::Cos(3*(Psi3WestEP[0])));
  2083. pCos3TPCvsRun[3][1]->Fill(event->runId(),TMath::Cos(4*(Psi3WestEP[0])));
  2084. pSin2TPCvsRun[0][1]->Fill(event->runId(),TMath::Sin(1*(Psi2WestEP[0])));
  2085. pSin2TPCvsRun[1][1]->Fill(event->runId(),TMath::Sin(2*(Psi2WestEP[0])));
  2086. pSin2TPCvsRun[2][1]->Fill(event->runId(),TMath::Sin(3*(Psi2WestEP[0])));
  2087. pSin2TPCvsRun[3][1]->Fill(event->runId(),TMath::Sin(4*(Psi2WestEP[0])));
  2088. pSin3TPCvsRun[0][1]->Fill(event->runId(),TMath::Sin(1*(Psi3WestEP[0])));
  2089. pSin3TPCvsRun[1][1]->Fill(event->runId(),TMath::Sin(2*(Psi3WestEP[0])));
  2090. pSin3TPCvsRun[2][1]->Fill(event->runId(),TMath::Sin(3*(Psi3WestEP[0])));
  2091. pSin3TPCvsRun[3][1]->Fill(event->runId(),TMath::Sin(4*(Psi3WestEP[0])));
  2092. }
  2093. if (fQcountEast[0] > cutNcountQvGap &&
  2094. fQcountWest[0] > cutNcountQvGap)
  2095. {
  2096. hTPCEastPsiShift[0]->Fill(GetEPAngle(Psi2EastEP[0],2), event->cent16());
  2097. hTPCWestPsiShift[0]->Fill(GetEPAngle(Psi2WestEP[0],2), event->cent16());
  2098. }
  2099. // 3rd harmonics TPC
  2100. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  2101. {
  2102. hTPCFullQxRaw[1]->Fill(fQv3FullEPRaw.X(),event->cent16());
  2103. hTPCFullQyRaw[1]->Fill(fQv3FullEPRaw.X(),event->cent16());
  2104. hTPCFullPsiRaw[1]->Fill(GetEPAngle(fQv3FullEPRaw,3), event->cent16());
  2105. }
  2106. if (fQcountEast[0] > cutNcountQvGap &&
  2107. fQcountWest[0] > cutNcountQvGap)
  2108. {
  2109. hTPCEastQxRaw[1]->Fill(fQv3EastEPRaw[0].X(),event->cent16());
  2110. hTPCEastQyRaw[1]->Fill(fQv3EastEPRaw[0].Y(),event->cent16());
  2111. hTPCWestQxRaw[1]->Fill(fQv3WestEPRaw[0].X(),event->cent16());
  2112. hTPCWestQyRaw[1]->Fill(fQv3WestEPRaw[0].Y(),event->cent16());
  2113. hTPCEastPsiRaw[1]->Fill(GetEPAngle(fQv3EastEPRaw[0],3), event->cent16());
  2114. hTPCWestPsiRaw[1]->Fill(GetEPAngle(fQv3WestEPRaw[0],3), event->cent16());
  2115. }
  2116. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  2117. {
  2118. hTPCFullQxRec[1]->Fill(fQv3FullEP.X(),event->cent16());
  2119. hTPCFullQyRec[1]->Fill(fQv3FullEP.X(),event->cent16());
  2120. hTPCFullPsiRec[1]->Fill(GetEPAngle(fQv3FullEP,3), event->cent16());
  2121. }
  2122. if (fQcountEast[0] > cutNcountQvGap &&
  2123. fQcountWest[0] > cutNcountQvGap)
  2124. {
  2125. hTPCEastQxRec[1]->Fill(fQv3EastEP[0].X(),event->cent16());
  2126. hTPCEastQyRec[1]->Fill(fQv3EastEP[0].Y(),event->cent16());
  2127. hTPCWestQxRec[1]->Fill(fQv3WestEP[0].X(),event->cent16());
  2128. hTPCWestQyRec[1]->Fill(fQv3WestEP[0].Y(),event->cent16());
  2129. hTPCEastPsiRec[1]->Fill(GetEPAngle(fQv3EastEP[0],3), event->cent16());
  2130. hTPCWestPsiRec[1]->Fill(GetEPAngle(fQv3WestEP[0],3), event->cent16());
  2131. }
  2132. if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
  2133. {
  2134. hTPCFullPsiShift[1]->Fill(GetEPAngle(Psi3FullEP,3), event->cent16());
  2135. }
  2136. if (fQcountEast[0] > cutNcountQvGap &&
  2137. fQcountWest[0] > cutNcountQvGap)
  2138. {
  2139. hTPCEastPsiShift[1]->Fill(GetEPAngle(Psi3EastEP[0],3), event->cent16());
  2140. hTPCWestPsiShift[1]->Fill(GetEPAngle(Psi3WestEP[0],3), event->cent16());
  2141. }
  2142. } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
  2143. femtoReader->Finish();
  2144. TFile *output = new TFile(outFile, "recreate");
  2145. for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
  2146. {
  2147. // Raw
  2148. hBBCFullQxRaw[iHarm]->Write();
  2149. hBBCFullQyRaw[iHarm]->Write();
  2150. hBBCEastQxRaw[iHarm]->Write();
  2151. hBBCEastQyRaw[iHarm]->Write();
  2152. hBBCWestQxRaw[iHarm]->Write();
  2153. hBBCWestQyRaw[iHarm]->Write();
  2154. hBBCFullPsiRaw[iHarm]->Write();
  2155. hBBCEastPsiRaw[iHarm]->Write();
  2156. hBBCWestPsiRaw[iHarm]->Write();
  2157. // w/ Gain correction
  2158. hBBCFullQxGain[iHarm]->Write();
  2159. hBBCFullQyGain[iHarm]->Write();
  2160. hBBCEastQxGain[iHarm]->Write();
  2161. hBBCEastQyGain[iHarm]->Write();
  2162. hBBCWestQxGain[iHarm]->Write();
  2163. hBBCWestQyGain[iHarm]->Write();
  2164. hBBCFullPsiGain[iHarm]->Write();
  2165. hBBCEastPsiGain[iHarm]->Write();
  2166. hBBCWestPsiGain[iHarm]->Write();
  2167. // w/ Gain + recentering
  2168. hBBCFullQxRec[iHarm]->Write();
  2169. hBBCFullQyRec[iHarm]->Write();
  2170. hBBCEastQxRec[iHarm]->Write();
  2171. hBBCEastQyRec[iHarm]->Write();
  2172. hBBCWestQxRec[iHarm]->Write();
  2173. hBBCWestQyRec[iHarm]->Write();
  2174. hBBCFullPsiRec[iHarm]->Write();
  2175. hBBCEastPsiRec[iHarm]->Write();
  2176. hBBCWestPsiRec[iHarm]->Write();
  2177. // w/ Gain + recentering + shift
  2178. hBBCFullPsiShift[iHarm]->Write();
  2179. hBBCEastPsiShift[iHarm]->Write();
  2180. hBBCWestPsiShift[iHarm]->Write();
  2181. }
  2182. for (int iHarm=0; iHarm<fNharmonics; iHarm++)
  2183. {
  2184. // Raw
  2185. hTPCFullQxRaw[iHarm]->Write();
  2186. hTPCFullQyRaw[iHarm]->Write();
  2187. hTPCEastQxRaw[iHarm]->Write();
  2188. hTPCEastQyRaw[iHarm]->Write();
  2189. hTPCWestQxRaw[iHarm]->Write();
  2190. hTPCWestQyRaw[iHarm]->Write();
  2191. hTPCFullPsiRaw[iHarm]->Write();
  2192. hTPCEastPsiRaw[iHarm]->Write();
  2193. hTPCWestPsiRaw[iHarm]->Write();
  2194. // w/ recentering
  2195. hTPCFullQxRec[iHarm]->Write();
  2196. hTPCFullQyRec[iHarm]->Write();
  2197. hTPCEastQxRec[iHarm]->Write();
  2198. hTPCEastQyRec[iHarm]->Write();
  2199. hTPCWestQxRec[iHarm]->Write();
  2200. hTPCWestQyRec[iHarm]->Write();
  2201. hTPCFullPsiRec[iHarm]->Write();
  2202. hTPCEastPsiRec[iHarm]->Write();
  2203. hTPCWestPsiRec[iHarm]->Write();
  2204. // w/ recentering + shift
  2205. hTPCFullPsiShift[iHarm]->Write();
  2206. hTPCEastPsiShift[iHarm]->Write();
  2207. hTPCWestPsiShift[iHarm]->Write();
  2208. }
  2209. for (int i = 0; i < 2; i++)
  2210. {
  2211. //Eventwise
  2212. hVtxXY[i]->Write();
  2213. hVtxZ[i]->Write();
  2214. hVtxVpdZ[i]->Write();
  2215. hVpdZ[i]->Write();
  2216. pVtxXvsRun[i]->Write();
  2217. pVtxYvsRun[i]->Write();
  2218. pVtxZvsRun[i]->Write();
  2219. pVpdZvsRun[i]->Write();
  2220. pVtxZVpdZvsRun[i]->Write();
  2221. pRefMultvsRun[i]->Write();
  2222. pTofMultvsRun[i]->Write();
  2223. pTofMatchvsRun[i]->Write();
  2224. pAdcvsRun[i]->Write();
  2225. pNtrPrimvsRun[i]->Write();
  2226. pNtrGlobvsRun[i]->Write();
  2227. //Trackwise
  2228. hNhitsFit[i]->Write();
  2229. hNhitsPoss[i]->Write();
  2230. hNhitsRatio[i]->Write();
  2231. hPt[i]->Write();
  2232. hPhi[i]->Write();
  2233. hPtot[i]->Write();
  2234. hEta[i]->Write();
  2235. hDCA[i]->Write();
  2236. pPtvsRun[i]->Write();
  2237. pEtavsRun[i]->Write();
  2238. pPhivsRun[i]->Write();
  2239. pNhitsFitvsRun[i]->Write();
  2240. pNhitsPossvsRun[i]->Write();
  2241. pDCAvsRun[i]->Write();
  2242. pdEdxvsRun[i]->Write();
  2243. pQx2TPCvsRun[i]->Write();
  2244. pQy2TPCvsRun[i]->Write();
  2245. pQx3TPCvsRun[i]->Write();
  2246. pQy3TPCvsRun[i]->Write();
  2247. pQx2BBCvsRun[i]->Write();
  2248. pQy2BBCvsRun[i]->Write();
  2249. for (int iCent=0;iCent<2*fNharmonics;iCent++)
  2250. {
  2251. pCos2TPCvsRun[iCent][i]->Write();
  2252. pSin2TPCvsRun[iCent][i]->Write();
  2253. pCos3TPCvsRun[iCent][i]->Write();
  2254. pSin3TPCvsRun[iCent][i]->Write();
  2255. pCosBBCvsRun[iCent][i]->Write();
  2256. pSinBBCvsRun[iCent][i]->Write();
  2257. }
  2258. //PID
  2259. hnSigM2_lowpt[i]->Write();
  2260. hnSigM2_highpt[i]->Write();
  2261. for (int iPID = 0; iPID < 4; iPID++)
  2262. {
  2263. hdEdxQp[i][iPID]->Write();
  2264. hM2Qp[i][iPID]->Write();
  2265. hNsigQp[i][iPID]->Write();
  2266. hM2pt[i][iPID]->Write();
  2267. hNsigpt[i][iPID]->Write();
  2268. }
  2269. }
  2270. for (int iPID = 0; iPID < 4; iPID++)
  2271. {
  2272. for (int iCent=0;iCent<16;iCent++)
  2273. {
  2274. h_pt_spectra[iPID][iCent]->Write();
  2275. }
  2276. }
  2277. hNevVsRun->Write();
  2278. output->Close();
  2279. std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
  2280. << std::endl;
  2281. timer.Stop();
  2282. timer.Print();
  2283. }
  2284. Bool_t isGoodEvent(StFemtoEvent *const &event)
  2285. {
  2286. if (!event)
  2287. return false;
  2288. if (event == nullptr)
  2289. return false;
  2290. if (event->primaryVertex().Perp() > cutVtxR)
  2291. return false;
  2292. if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy))
  2293. return false;
  2294. if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz)
  2295. return false;
  2296. return true;
  2297. }
  2298. Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
  2299. {
  2300. if (!track)
  2301. return false;
  2302. // if (!track->isPrimary()) return false;
  2303. if (track->nHitsFit() < cutNhits)
  2304. return false;
  2305. if (track->nHitsPoss() <= cutNhitsPoss)
  2306. return false;
  2307. if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
  2308. return false;
  2309. if (TMath::Abs(track->eta()) >= cutEta)
  2310. return false;
  2311. if (track->pt() <= cutPtMin.at(_energy))
  2312. return false;
  2313. if (track->pt() > cutPtMax)
  2314. return false;
  2315. if (track->ptot() > cutPMax)
  2316. return false;
  2317. if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
  2318. return false;
  2319. return true;
  2320. }
  2321. Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
  2322. {
  2323. if (!track)
  2324. return false;
  2325. // if (!track->isPrimary()) return false;
  2326. if (track->nHitsFit() < cutNhits)
  2327. return false;
  2328. if (track->nHitsPoss() <= cutNhitsPoss)
  2329. return false;
  2330. if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
  2331. return false;
  2332. if (TMath::Abs(track->eta()) >= cutEta)
  2333. return false;
  2334. if (track->pt() <= cutPtMin.at(_energy))
  2335. return false;
  2336. //if (track->pt() > cutPtMax) return false;
  2337. if (track->ptot() > cutPMax)
  2338. return false;
  2339. if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
  2340. return false;
  2341. return true;
  2342. }
  2343. Double_t GetWeight(StFemtoTrack *const &track)
  2344. {
  2345. Double_t weight;
  2346. if (track->pt() <= cutPtWeightEP)
  2347. {
  2348. weight = track->pt();
  2349. }
  2350. else
  2351. {
  2352. weight = cutPtWeightEP;
  2353. }
  2354. return weight;
  2355. }
  2356. TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
  2357. {
  2358. TVector2 qv(0., 0.);
  2359. qv.Set(TMath::Cos(_harm * track->phi()), TMath::Sin(_harm * track->phi()));
  2360. return qv;
  2361. }
  2362. Double_t AngleShift(Double_t Psi, Double_t order)
  2363. {
  2364. Double_t PsiCorr = Psi;
  2365. if (PsiCorr > TMath::Pi() / order)
  2366. {
  2367. PsiCorr = Psi - 2. * TMath::Pi() / order;
  2368. }
  2369. if (PsiCorr < -1. * TMath::Pi() / order)
  2370. {
  2371. PsiCorr = Psi + 2. * TMath::Pi() / order;
  2372. }
  2373. return PsiCorr;
  2374. }
  2375. Int_t GetVzBin(Double_t vtxZ)
  2376. {
  2377. for (Int_t i = 0; i < fArraySize; i++)
  2378. {
  2379. if (vtxZ >= VtxBins[i] && vtxZ < VtxBins[i + 1])
  2380. return i;
  2381. }
  2382. return -1;
  2383. }
  2384. //--------------------------------------------------------------------------------------------------------------------------//
  2385. //this function simply connects the gain values read in to the BBC azimuthal distribution
  2386. //since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
  2387. //I prefer assigning the angle between the tiles thus "greying out" the adcs.
  2388. //Others have assigned all of the adc to one (exclusive) or the the other.
  2389. Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId)
  2390. {
  2391. //float GetPhiInBBC(int eastWest, int bbcN) { //tileId=0 to 23
  2392. const float Pi = TMath::Pi();
  2393. const float Phi_div = Pi / 6;
  2394. float bbc_phi = Phi_div;
  2395. switch(tileId) {
  2396. case 0: bbc_phi = 3.*Phi_div;
  2397. break;
  2398. case 1: bbc_phi = Phi_div;
  2399. break;
  2400. case 2: bbc_phi = -1.*Phi_div;
  2401. break;
  2402. case 3: bbc_phi = -3.*Phi_div;
  2403. break;
  2404. case 4: bbc_phi = -5.*Phi_div;
  2405. break;
  2406. case 5: bbc_phi = 5.*Phi_div;
  2407. break;
  2408. //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
  2409. case 6: bbc_phi = 3.*Phi_div;
  2410. break;
  2411. case 7: bbc_phi = 3.*Phi_div;
  2412. break;
  2413. case 8: bbc_phi = Phi_div;
  2414. break;
  2415. case 9: bbc_phi = 0.;
  2416. break;
  2417. case 10: bbc_phi = -1.*Phi_div;
  2418. break;
  2419. //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div; //tiles 13 and 15 are gained together
  2420. case 11: bbc_phi = -3.*Phi_div;
  2421. break;
  2422. case 12: bbc_phi = -3.*Phi_div;
  2423. break;
  2424. case 13: bbc_phi = -5.*Phi_div;
  2425. break;
  2426. case 14: bbc_phi = Pi;
  2427. break;
  2428. case 15: bbc_phi = 5.*Phi_div;
  2429. break;
  2430. }
  2431. //if we're looking at the east BBC we need to flip around x in the STAR coordinates,
  2432. //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the
  2433. if(0 == eastWest){
  2434. if (bbc_phi > -0.001){ //this is not a >= since we are talking about finite adcs -- not to important
  2435. bbc_phi = Pi - bbc_phi;
  2436. }
  2437. else {
  2438. bbc_phi= -Pi - bbc_phi;
  2439. }
  2440. }
  2441. if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
  2442. if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;
  2443. return bbc_phi;
  2444. }
  2445. Double_t GetEPAngle(TVector2 qvect, Double_t order)
  2446. {
  2447. Double_t angle = -9999.;
  2448. if(qvect.Mod())
  2449. {
  2450. angle = qvect.Phi();
  2451. if(angle < 0.0) {angle += 2.*TMath::Pi()/(order);}
  2452. if(angle > 2.*TMath::Pi()/(order)) {angle -= 2.*TMath::Pi()/(order);}
  2453. }
  2454. return angle;
  2455. }
  2456. Double_t GetEPAngle(Double_t angle, Double_t order)
  2457. {
  2458. if(angle < 0.0) {angle += 2.*TMath::Pi()/(order);}
  2459. if(angle > 2.*TMath::Pi()/(order)) {angle -= 2.*TMath::Pi()/(order);}
  2460. return angle;
  2461. }
  2462. Bool_t isGoodPID(StFemtoTrack *const &track)
  2463. {
  2464. if (!track->isTofTrack()) return false;
  2465. if (track->massSqr() < cutMass2Min) return false;
  2466. //NhitsDedx cut applied in StFemtoDst?
  2467. // ToFYLocal cut applied in StFemtoDst?
  2468. return true;
  2469. }