12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406140714081409141014111412141314141415141614171418141914201421142214231424142514261427142814291430143114321433143414351436143714381439144014411442144314441445144614471448144914501451145214531454145514561457145814591460146114621463146414651466146714681469147014711472147314741475147614771478147914801481148214831484148514861487148814891490149114921493149414951496149714981499150015011502150315041505150615071508150915101511151215131514151515161517151815191520152115221523152415251526152715281529153015311532153315341535153615371538153915401541154215431544154515461547154815491550155115521553155415551556155715581559156015611562156315641565156615671568156915701571157215731574157515761577157815791580158115821583158415851586158715881589159015911592159315941595159615971598159916001601160216031604160516061607160816091610161116121613161416151616161716181619162016211622162316241625162616271628162916301631163216331634163516361637163816391640164116421643164416451646164716481649165016511652165316541655165616571658165916601661166216631664166516661667166816691670167116721673167416751676167716781679168016811682168316841685168616871688168916901691169216931694169516961697169816991700170117021703170417051706170717081709171017111712171317141715171617171718171917201721172217231724172517261727172817291730173117321733173417351736173717381739174017411742174317441745174617471748174917501751175217531754175517561757175817591760176117621763176417651766176717681769177017711772177317741775177617771778177917801781178217831784178517861787178817891790179117921793179417951796179717981799180018011802180318041805180618071808180918101811181218131814181518161817181818191820182118221823182418251826182718281829183018311832183318341835183618371838183918401841184218431844184518461847184818491850185118521853185418551856185718581859186018611862186318641865186618671868186918701871187218731874187518761877187818791880188118821883188418851886188718881889189018911892189318941895189618971898189919001901190219031904190519061907190819091910191119121913191419151916191719181919192019211922192319241925192619271928192919301931193219331934193519361937193819391940194119421943194419451946194719481949195019511952195319541955195619571958195919601961196219631964196519661967196819691970197119721973197419751976197719781979198019811982198319841985198619871988198919901991199219931994199519961997199819992000200120022003200420052006200720082009201020112012201320142015201620172018201920202021202220232024202520262027202820292030203120322033203420352036203720382039204020412042204320442045204620472048204920502051205220532054205520562057205820592060206120622063206420652066206720682069207020712072207320742075207620772078207920802081208220832084208520862087208820892090209120922093209420952096209720982099210021012102210321042105210621072108210921102111211221132114211521162117211821192120212121222123212421252126212721282129213021312132213321342135213621372138213921402141214221432144214521462147214821492150215121522153215421552156215721582159216021612162216321642165216621672168216921702171217221732174217521762177217821792180218121822183218421852186218721882189219021912192219321942195219621972198219922002201220222032204220522062207220822092210221122122213221422152216221722182219222022212222222322242225222622272228222922302231223222332234223522362237223822392240224122422243224422452246224722482249225022512252225322542255225622572258225922602261226222632264226522662267226822692270227122722273227422752276227722782279228022812282228322842285228622872288228922902291229222932294229522962297229822992300230123022303230423052306230723082309231023112312231323142315231623172318231923202321232223232324232523262327232823292330233123322333233423352336233723382339234023412342234323442345234623472348234923502351235223532354235523562357235823592360236123622363236423652366236723682369237023712372237323742375237623772378237923802381238223832384238523862387238823892390239123922393239423952396239723982399240024012402240324042405240624072408240924102411241224132414241524162417241824192420242124222423242424252426242724282429243024312432243324342435243624372438243924402441244224432444244524462447244824492450245124522453245424552456245724582459246024612462246324642465246624672468246924702471247224732474247524762477247824792480248124822483248424852486248724882489249024912492249324942495249624972498249925002501250225032504250525062507250825092510251125122513251425152516251725182519252025212522252325242525252625272528252925302531253225332534253525362537253825392540254125422543254425452546254725482549255025512552255325542555255625572558255925602561256225632564256525662567256825692570257125722573257425752576257725782579258025812582258325842585258625872588258925902591259225932594259525962597259825992600260126022603260426052606260726082609261026112612261326142615261626172618261926202621262226232624262526262627262826292630263126322633263426352636263726382639264026412642264326442645264626472648264926502651265226532654265526562657265826592660266126622663266426652666266726682669267026712672267326742675267626772678267926802681268226832684268526862687268826892690269126922693269426952696269726982699270027012702270327042705270627072708270927102711271227132714271527162717 |
- /**
- * \brief Example of how to read a file (list of files) using StFemtoEvent classes
- *
- * RunFemtoDstAnalyzer.C is an example of reading FemtoDst format.
- * One can use either FemtoDst file or a list of femtoDst files (inFile.lis or
- * inFile.list) as an input, and preform physics analysis
- *
- * \author Grigory Nigmatkulov
- * \date May 29, 2018
- */
- // This is needed for calling standalone classes
- #define _VANILLA_ROOT_
- // C++ headers
- #include <iostream>
- #include <vector>
- #include <map>
- // ROOT headers
- #include "TROOT.h"
- #include "TFile.h"
- #include "TChain.h"
- #include "TVector2.h"
- #include "TVector3.h"
- #include "TTree.h"
- #include "TSystem.h"
- #include "TH1.h"
- #include "TH2.h"
- #include "TMath.h"
- #include "TProfile2D.h"
- #include "TStopwatch.h"
- #include "TRandom3.h"
- #include "TVector2.h"
- // FemtoDst headers
- #include "StFemtoDstReader.h"
- #include "StFemtoDst.h"
- #include "StFemtoEvent.h"
- #include "StFemtoTrack.h"
- // Load libraries (for ROOT_VERSTION_CODE >= 393215)
- #if ROOT_VERSION_CODE >= ROOT_VERSION(6, 0, 0)
- R__LOAD_LIBRARY(/ mnt / pool / rhic / 4 / parfenovpeter / STAR / build / libStFemtoDst.so)
- #endif
- #include "Constants.h"
- // inFile - is a name of name.FemtoDst.root file or a name
- // of a name.lis(t) files, that contains a list of
- // name1.FemtoDst.root, name2.FemtoDst.root, ... files
- //Used functions (see them below)
- Bool_t isGoodEvent(StFemtoEvent *const &event);
- Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
- Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx);
- Double_t GetWeight(StFemtoTrack *const &track);
- TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm);
- Double_t AngleShift(Double_t Psi, Double_t order);
- Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);
- Int_t GetVzBin(Double_t vtxZ);
- Double_t GetEPAngle(TVector2 qvect, Double_t order);
- Double_t GetEPAngle(Double_t angle, Double_t order);
- Bool_t isGoodPID(StFemtoTrack *const &track);
- //_________________
- void FemtoDstAnalyzer_QA(const Char_t *inFile = "st_physics_12150008_raw_4030001.femtoDst.root",
- const Char_t *outFile = "./oShiftTest.root",
- const Char_t *reCentFile = "",
- const Char_t *shiftFile = "",
- const Char_t *gainBBC = "",
- const Char_t *reCentFileBBC = "",
- const Char_t *shiftFileBBC = "")
- {
- std::cout << "Hi! Lets do some physics, Master!" << std::endl;
- TStopwatch timer;
- timer.Start();
- #if ROOT_VERSION_CODE < ROOT_VERSION(6, 0, 0)
- gSystem->Load("../libStFemtoDst.so");
- #endif
- StFemtoDstReader *femtoReader = new StFemtoDstReader(inFile);
- femtoReader->Init();
- // This is a way if you want to spead up IO
- std::cout << "Explicit read status for some branches" << std::endl;
- femtoReader->SetStatus("*", 0);
- femtoReader->SetStatus("Event", 1);
- femtoReader->SetStatus("Track", 1);
- std::cout << "Status has been set" << std::endl;
- std::cout << "Now I know what to read, Master!" << std::endl;
- if (!femtoReader->chain())
- {
- std::cout << "No chain has been found." << std::endl;
- }
- Long64_t eventsInTree = femtoReader->tree()->GetEntries();
- std::cout << "eventsInTree: " << eventsInTree << std::endl;
- Long64_t events2read = femtoReader->chain()->GetEntries();
- std::cout << "Number of events to read: " << events2read << std::endl;
- //Read ReCentering <Q> vectors
- TFile *fiReCent = new TFile(reCentFile, "read");
- std::map<Double_t, Double_t> p_Qx_FULL_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_Qy_FULL_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_Qx_FULL_SP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_Qy_FULL_SP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_Qx_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Qy_EAST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Qx_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Qy_EAST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Qx_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Qy_WEST_EP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Qx_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Qy_WEST_SP[fNharmonics][fArraySize][fNCent][NEtaGaps];
- TProfile2D *prof;
- double binCont;
- std::cout << "Reading file for TPC recentering..." << std::endl;
- for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
- {
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_EP%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qx_FULL_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_EP%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qy_FULL_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_FULL_SP%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qx_FULL_SP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_FULL_SP%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qy_FULL_SP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- for (int iEtaGap = 0; iEtaGap < NEtaGaps; iEtaGap++)
- {
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qx_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qy_EAST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qx_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_EAST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qy_EAST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qx_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_EP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qy_WEST_EP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%ix_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qx_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCent->Get(Form("p_Q%iy_WEST_SP%i_gap%i", iHarm, iVtx, iEtaGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Qy_WEST_SP[iHarm][iVtx][iCent][iEtaGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- }
- }
- }
- //Read Shift coefficients
- TFile *fiShift = new TFile(shiftFile, "read");
- std::map<Double_t, Double_t> p_Cos2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Sin2_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Cos2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Sin2_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Cos2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin2_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Cos2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin2_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Cos2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin2_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Cos2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin2_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Cos3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Sin3_FULL_EP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Cos3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Sin3_FULL_SP[fArraySize][NShiftOrderMax][fNCent];
- std::map<Double_t, Double_t> p_Cos3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin3_EAST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Cos3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin3_EAST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Cos3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin3_WEST_EP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Cos3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::map<Double_t, Double_t> p_Sin3_WEST_SP[fArraySize][NShiftOrderMax][fNCent][NEtaGaps];
- std::cout << "Reading file for TPC flattening..." << std::endl;
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- for (int iShift = 0; iShift < NShiftOrderMax; iShift++)
- {
- std::cout << "iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos2_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 2, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin2_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos2_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 2, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin2_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos3_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_EP%i_shiftOrder_%i", 3, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin3_FULL_EP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos3_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_FULL_SP%i_shiftOrder_%i", 3, iVtx, iShift));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin3_FULL_SP[iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- for (int iGap = 0; iGap < NEtaGaps; iGap++)
- {
- // for 2nd harmonics
- // EAST
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin2_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin2_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- // WEST
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin2_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 2, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin2_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- // for 3rd harmonics
- // EAST
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin3_EAST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_EAST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin3_EAST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- // WEST
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_EP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin3_WEST_EP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Cos%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Cos3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShift->Get(Form("p_Sin%i_WEST_SP%i_shiftOrder_%i_gap%i", 3, iVtx, iShift, iGap));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_Sin3_WEST_SP[iVtx][iShift][iCent][iGap][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- }
- }
- }
- //Manage gain correction
- TProfile2D *BBCgain_East[fArraySize];
- TProfile2D *BBCgain_West[fArraySize];
- TFile *fiGain = new TFile(gainBBC, "read");
- fiGain->cd();
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- BBCgain_East[iVtx] = (TProfile2D *)fiGain->Get(Form("p_gainBBC_East_Vtx%i", iVtx));
- BBCgain_West[iVtx] = (TProfile2D *)fiGain->Get(Form("p_gainBBC_West_Vtx%i", iVtx));
- }
- //Manage recentering BBC
- std::map<Double_t, Double_t> p_BBC_Qx_Full_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_BBC_Qy_Full_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_BBC_Qx_East_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_BBC_Qy_East_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_BBC_Qx_West_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_BBC_Qy_West_EP[fNharmonics][fArraySize][fNCent];
- TFile *fiReCentBBC = new TFile(reCentFileBBC, "read");
- std::cout << "Reading file for BBC recentering..." << std::endl;
- for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
- {
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << std::endl;
- prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_Full_EP%i_vtx%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Qx_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_Full_EP%i_vtx%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Qy_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_East_EP%i_vtx%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Qx_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_East_EP%i_vtx%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Qy_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qx_West_EP%i_vtx%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Qx_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentBBC->Get(Form("p_BBC_Qy_West_EP%i_vtx%i", iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Qy_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- }
- }
- //Manage ShiftCorr BBC
- std::map<Double_t, Double_t> p_BBC_Cos_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_BBC_Sin_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_BBC_Cos_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_BBC_Sin_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_BBC_Cos_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_BBC_Sin_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- TFile *fiShiftBBC = new TFile(shiftFileBBC, "read");
- fiShiftBBC->cd();
- std::cout << "Reading file for BBC shift corrections..." << std::endl;
- for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- {
- for (int iVtx = 0; iVtx < fArraySize; iVtx++)
- {
- for (int iShift = 0; iShift < NShiftOrderMax * 3; iShift++)
- {
- std::cout << "iHarm: " << iHarm << " iVtx: " << iVtx << " iShiftOrder: " << iShift << std::endl;
- prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Cos_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_Full_EP%i_vtx%i", iShift, iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Sin_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Cos_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_East_EP%i_vtx%i", iShift, iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Sin_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Cos%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Cos_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftBBC->Get(Form("p_BBC_Sin%i_West_EP%i_vtx%i", iShift, iHarm, iVtx));
- for (int i = 0; i < fNBins; i++)
- {
- for (int iCent = 0; iCent < fNCent; iCent++)
- {
- binCont = prof->GetBinContent(prof->FindBin((Double_t)(runId_min + i), (Double_t)iCent));
- if (binCont == 0.)
- continue;
- p_BBC_Sin_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- }
- }
- }
- //Eventwise parameters
- //hHisto[0] - before event selection
- //hHisto[1] - after event selection
- TH2D *hVtxXY[2];
- TH1D *hVtxZ[2];
- TH1D *hVpdZ[2];
- TH1D *hVtxVpdZ[2];
- TProfile *pVtxZvsRun[2];
- TProfile *pVtxXvsRun[2];
- TProfile *pVtxYvsRun[2];
- TProfile *pVpdZvsRun[2];
- TProfile *pVtxZVpdZvsRun[2];
- TProfile *pRefMultvsRun[2];
- TProfile *pAdcvsRun[2];
- TProfile *pTofMultvsRun[2];
- TProfile *pTofMatchvsRun[2];
- TH2D *hNevVsRun = new TH2D(Form("hNevVsRun"),Form("hNevVsRun"),fNBins, runId_min - 0.5, runId_max + 0.5,16,-0.5,15.5);
- //BBC EP related parameters
- // Raw
- TH2D* hBBCFullQxRaw[fNharmonics+1];
- TH2D* hBBCFullQyRaw[fNharmonics+1];
- TH2D* hBBCEastQxRaw[fNharmonics+1];
- TH2D* hBBCEastQyRaw[fNharmonics+1];
- TH2D* hBBCWestQxRaw[fNharmonics+1];
- TH2D* hBBCWestQyRaw[fNharmonics+1];
- TH2D* hBBCFullPsiRaw[fNharmonics+1];
- TH2D* hBBCEastPsiRaw[fNharmonics+1];
- TH2D* hBBCWestPsiRaw[fNharmonics+1];
- // w/ Gain correction
- TH2D* hBBCFullQxGain[fNharmonics+1];
- TH2D* hBBCFullQyGain[fNharmonics+1];
- TH2D* hBBCEastQxGain[fNharmonics+1];
- TH2D* hBBCEastQyGain[fNharmonics+1];
- TH2D* hBBCWestQxGain[fNharmonics+1];
- TH2D* hBBCWestQyGain[fNharmonics+1];
- TH2D* hBBCFullPsiGain[fNharmonics+1];
- TH2D* hBBCEastPsiGain[fNharmonics+1];
- TH2D* hBBCWestPsiGain[fNharmonics+1];
- // w/ Gain + recentering
- TH2D* hBBCFullQxRec[fNharmonics+1];
- TH2D* hBBCFullQyRec[fNharmonics+1];
- TH2D* hBBCEastQxRec[fNharmonics+1];
- TH2D* hBBCEastQyRec[fNharmonics+1];
- TH2D* hBBCWestQxRec[fNharmonics+1];
- TH2D* hBBCWestQyRec[fNharmonics+1];
- TH2D* hBBCFullPsiRec[fNharmonics+1];
- TH2D* hBBCEastPsiRec[fNharmonics+1];
- TH2D* hBBCWestPsiRec[fNharmonics+1];
- // w/ Gain + recentering + shift
- TH2D* hBBCFullPsiShift[fNharmonics+1];
- TH2D* hBBCEastPsiShift[fNharmonics+1];
- TH2D* hBBCWestPsiShift[fNharmonics+1];
- //ZDC EP related parameters
- // Raw
- TH2D* hZDCFullQxRaw;
- TH2D* hZDCFullQyRaw;
- TH2D* hZDCEastQxRaw;
- TH2D* hZDCEastQyRaw;
- TH2D* hZDCWestQxRaw;
- TH2D* hZDCWestQyRaw;
- TH2D* hZDCFullPsiRaw;
- TH2D* hZDCEastPsiRaw;
- TH2D* hZDCWestPsiRaw;
- // w/ recentering
- TH2D* hZDCFullQxRec;
- TH2D* hZDCFullQyRec;
- TH2D* hZDCEastQxRec;
- TH2D* hZDCEastQyRec;
- TH2D* hZDCWestQxRec;
- TH2D* hZDCWestQyRec;
- TH2D* hZDCFullPsiRec;
- TH2D* hZDCEastPsiRec;
- TH2D* hZDCWestPsiRec;
- // w/ recentering + shift
- TH2D* hZDCFullPsiShift;
- TH2D* hZDCEastPsiShift;
- TH2D* hZDCWestPsiShift;
- //TPC EP related parameters
- // Raw
- TH2D* hTPCFullQxRaw[fNharmonics];
- TH2D* hTPCFullQyRaw[fNharmonics];
- TH2D* hTPCEastQxRaw[fNharmonics];
- TH2D* hTPCEastQyRaw[fNharmonics];
- TH2D* hTPCWestQxRaw[fNharmonics];
- TH2D* hTPCWestQyRaw[fNharmonics];
- TH2D* hTPCFullPsiRaw[fNharmonics];
- TH2D* hTPCEastPsiRaw[fNharmonics];
- TH2D* hTPCWestPsiRaw[fNharmonics];
- // w/ recentering
- TH2D* hTPCFullQxRec[fNharmonics];
- TH2D* hTPCFullQyRec[fNharmonics];
- TH2D* hTPCEastQxRec[fNharmonics];
- TH2D* hTPCEastQyRec[fNharmonics];
- TH2D* hTPCWestQxRec[fNharmonics];
- TH2D* hTPCWestQyRec[fNharmonics];
- TH2D* hTPCFullPsiRec[fNharmonics];
- TH2D* hTPCEastPsiRec[fNharmonics];
- TH2D* hTPCWestPsiRec[fNharmonics];
- // w/ recentering + shift
- TH2D* hTPCFullPsiShift[fNharmonics];
- TH2D* hTPCEastPsiShift[fNharmonics];
- TH2D* hTPCWestPsiShift[fNharmonics];
- TProfile *pQx2TPCvsRun[2]; // 0 - East, 1 - West
- TProfile *pQy2TPCvsRun[2]; // 0 - East, 1 - West
- TProfile *pQx3TPCvsRun[2]; // 0 - East, 1 - West
- TProfile *pQy3TPCvsRun[2]; // 0 - East, 1 - West
- TProfile *pQx2BBCvsRun[2]; // 0 - East, 1 - West
- TProfile *pQy2BBCvsRun[2]; // 0 - East, 1 - West
- TProfile *pQx1ZDCvsRun[2]; // 0 - East, 1 - West
- TProfile *pQy1ZDCvsRun[2]; // 0 - East, 1 - West
- TProfile *pCos2TPCvsRun[2*fNharmonics][2];
- TProfile *pSin2TPCvsRun[2*fNharmonics][2];
- TProfile *pCos3TPCvsRun[2*fNharmonics][2];
- TProfile *pSin3TPCvsRun[2*fNharmonics][2];
- TProfile *pCosBBCvsRun[2*fNharmonics][2];
- TProfile *pSinBBCvsRun[2*fNharmonics][2];
- TProfile *pCosZDCvsRun[2];
- TProfile *pSinZDCvsRun[2];
- //Trackwise parameters
- //hHisto[0] - before track selection (but after event selection)
- //hHisto[1] - after track selection
- TH1D *hNhitsFit[2];
- TH1D *hNhitsPoss[2];
- TH1D *hNhitsRatio[2];
- TH1D *hPt[2];
- TH1D *hPtot[2];
- TH1D *hPhi[2];
- TH1D *hEta[2];
- TH1D *hDCA[2];
- TProfile *pPtvsRun[2];
- TProfile *pEtavsRun[2];
- TProfile *pPhivsRun[2];
- TProfile *pNhitsFitvsRun[2];
- TProfile *pNhitsPossvsRun[2];
- TProfile *pDCAvsRun[2];
- TProfile *pdEdxvsRun[2];
- TProfile *pNtrPrimvsRun[2];
- TProfile *pNtrGlobvsRun[2];
- //PID related parameters
- //hHisto[0] - before any PID (but after track selection)
- //hHisto[1] - after PID
- //second [i] is for particle species: 0-all, 1-pion, 2-kaon, 3-proton
- TH2D *hnSigM2_lowpt[2];
- TH2D *hnSigM2_highpt[2];
- TH2D *hdEdxQp[2][4];
- TH2D *hM2Qp[2][4];
- TH2D *hNsigQp[2][4];
- TH2D *hM2pt[2][4];
- TH2D *hNsigpt[2][4];
- TProfile2D *BBCgain[2][24];
- TH2D *h_pt_spectra[4][16]; //for <pT> calc 0 - pions, 1 - kaons, 2 - protons, 3 - charged hadrons
- //Initialization
- //ZDC
- // Raw
- hZDCFullQxRaw = new TH2D(Form("hZDCFullQxRaw%i",0), Form("hZDCFullQxRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCFullQyRaw = new TH2D(Form("hZDCFullQyRaw%i",0), Form("hZDCFullQyRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCEastQxRaw = new TH2D(Form("hZDCEastQxRaw%i",0), Form("hZDCEastQxRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCEastQyRaw = new TH2D(Form("hZDCEastQyRaw%i",0), Form("hZDCEastQyRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCWestQxRaw = new TH2D(Form("hZDCWestQxRaw%i",0), Form("hZDCWestQxRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCWestQyRaw = new TH2D(Form("hZDCWestQyRaw%i",0), Form("hZDCWestQyRaw%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCFullPsiRaw = new TH2D(Form("hZDCFullPsiRaw%i",0), Form("hZDCFullPsiRaw%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- hZDCEastPsiRaw = new TH2D(Form("hZDCEastPsiRaw%i",0), Form("hZDCEastPsiRaw%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- hZDCWestPsiRaw = new TH2D(Form("hZDCWestPsiRaw%i",0), Form("hZDCWestPsiRaw%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- // Rec
- hZDCFullQxRec = new TH2D(Form("hZDCFullQxRec%i",0), Form("hZDCFullQxRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCFullQyRec = new TH2D(Form("hZDCFullQyRec%i",0), Form("hZDCFullQyRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCEastQxRec = new TH2D(Form("hZDCEastQxRec%i",0), Form("hZDCEastQxRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCEastQyRec = new TH2D(Form("hZDCEastQyRec%i",0), Form("hZDCEastQyRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCWestQxRec = new TH2D(Form("hZDCWestQxRec%i",0), Form("hZDCWestQxRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCWestQyRec = new TH2D(Form("hZDCWestQyRec%i",0), Form("hZDCWestQyRec%i",0), 200, -1.2,1.2, 16, -0.5, 15.5);
- hZDCFullPsiRec = new TH2D(Form("hZDCFullPsiRec%i",0), Form("hZDCFullPsiRec%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- hZDCEastPsiRec = new TH2D(Form("hZDCEastPsiRec%i",0), Form("hZDCEastPsiRec%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- hZDCWestPsiRec = new TH2D(Form("hZDCWestPsiRec%i",0), Form("hZDCWestPsiRec%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- // Shift
- hZDCFullPsiShift = new TH2D(Form("hZDCFullPsiShift%i",0), Form("hZDCFullPsiShift%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- hZDCEastPsiShift = new TH2D(Form("hZDCEastPsiShift%i",0), Form("hZDCEastPsiShift%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
- hZDCWestPsiShift = new TH2D(Form("hZDCWestPsiShift%i",0), Form("hZDCWestPsiShift%i",0), 720, 0., 4.*TMath::Pi()/(0+1.), 16, -0.5, 15.5);
-
- //BBC
- for (int iHarm = 0; iHarm < fNharmonics+1; iHarm++)
- {
- // Raw
- hBBCFullQxRaw[iHarm] = new TH2D(Form("hBBCFullQxRaw%i",iHarm), Form("hBBCFullQxRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCFullQyRaw[iHarm] = new TH2D(Form("hBBCFullQyRaw%i",iHarm), Form("hBBCFullQyRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCEastQxRaw[iHarm] = new TH2D(Form("hBBCEastQxRaw%i",iHarm), Form("hBBCEastQxRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCEastQyRaw[iHarm] = new TH2D(Form("hBBCEastQyRaw%i",iHarm), Form("hBBCEastQyRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCWestQxRaw[iHarm] = new TH2D(Form("hBBCWestQxRaw%i",iHarm), Form("hBBCWestQxRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCWestQyRaw[iHarm] = new TH2D(Form("hBBCWestQyRaw%i",iHarm), Form("hBBCWestQyRaw%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCFullPsiRaw[iHarm] = new TH2D(Form("hBBCFullPsiRaw%i",iHarm), Form("hBBCFullPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCEastPsiRaw[iHarm] = new TH2D(Form("hBBCEastPsiRaw%i",iHarm), Form("hBBCEastPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCWestPsiRaw[iHarm] = new TH2D(Form("hBBCWestPsiRaw%i",iHarm), Form("hBBCWestPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- // Gain
- hBBCFullQxGain[iHarm] = new TH2D(Form("hBBCFullQxGain%i",iHarm), Form("hBBCFullQxGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCFullQyGain[iHarm] = new TH2D(Form("hBBCFullQyGain%i",iHarm), Form("hBBCFullQyGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCEastQxGain[iHarm] = new TH2D(Form("hBBCEastQxGain%i",iHarm), Form("hBBCEastQxGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCEastQyGain[iHarm] = new TH2D(Form("hBBCEastQyGain%i",iHarm), Form("hBBCEastQyGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCWestQxGain[iHarm] = new TH2D(Form("hBBCWestQxGain%i",iHarm), Form("hBBCWestQxGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCWestQyGain[iHarm] = new TH2D(Form("hBBCWestQyGain%i",iHarm), Form("hBBCWestQyGain%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCFullPsiGain[iHarm] = new TH2D(Form("hBBCFullPsiGain%i",iHarm), Form("hBBCFullPsiGain%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCEastPsiGain[iHarm] = new TH2D(Form("hBBCEastPsiGain%i",iHarm), Form("hBBCEastPsiGain%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCWestPsiGain[iHarm] = new TH2D(Form("hBBCWestPsiGain%i",iHarm), Form("hBBCWestPsiGain%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- // Rec
- hBBCFullQxRec[iHarm] = new TH2D(Form("hBBCFullQxRec%i",iHarm), Form("hBBCFullQxRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCFullQyRec[iHarm] = new TH2D(Form("hBBCFullQyRec%i",iHarm), Form("hBBCFullQyRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCEastQxRec[iHarm] = new TH2D(Form("hBBCEastQxRec%i",iHarm), Form("hBBCEastQxRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCEastQyRec[iHarm] = new TH2D(Form("hBBCEastQyRec%i",iHarm), Form("hBBCEastQyRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCWestQxRec[iHarm] = new TH2D(Form("hBBCWestQxRec%i",iHarm), Form("hBBCWestQxRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCWestQyRec[iHarm] = new TH2D(Form("hBBCWestQyRec%i",iHarm), Form("hBBCWestQyRec%i",iHarm), 200, -1.2,1.2, 16, -0.5, 15.5);
- hBBCFullPsiRec[iHarm] = new TH2D(Form("hBBCFullPsiRec%i",iHarm), Form("hBBCFullPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCEastPsiRec[iHarm] = new TH2D(Form("hBBCEastPsiRec%i",iHarm), Form("hBBCEastPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCWestPsiRec[iHarm] = new TH2D(Form("hBBCWestPsiRec%i",iHarm), Form("hBBCWestPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- // Shift
- hBBCFullPsiShift[iHarm] = new TH2D(Form("hBBCFullPsiShift%i",iHarm), Form("hBBCFullPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCEastPsiShift[iHarm] = new TH2D(Form("hBBCEastPsiShift%i",iHarm), Form("hBBCEastPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hBBCWestPsiShift[iHarm] = new TH2D(Form("hBBCWestPsiShift%i",iHarm), Form("hBBCWestPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- }
- //TPC
- for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
- {
- // Raw
- hTPCFullQxRaw[iHarm] = new TH2D(Form("hTPCFullQxRaw%i",iHarm), Form("hTPCFullQxRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCFullQyRaw[iHarm] = new TH2D(Form("hTPCFullQyRaw%i",iHarm), Form("hTPCFullQyRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCEastQxRaw[iHarm] = new TH2D(Form("hTPCEastQxRaw%i",iHarm), Form("hTPCEastQxRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCEastQyRaw[iHarm] = new TH2D(Form("hTPCEastQyRaw%i",iHarm), Form("hTPCEastQyRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCWestQxRaw[iHarm] = new TH2D(Form("hTPCWestQxRaw%i",iHarm), Form("hTPCWestQxRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCWestQyRaw[iHarm] = new TH2D(Form("hTPCWestQyRaw%i",iHarm), Form("hTPCWestQyRaw%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCFullPsiRaw[iHarm] = new TH2D(Form("hTPCFullPsiRaw%i",iHarm), Form("hTPCFullPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hTPCEastPsiRaw[iHarm] = new TH2D(Form("hTPCEastPsiRaw%i",iHarm), Form("hTPCEastPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hTPCWestPsiRaw[iHarm] = new TH2D(Form("hTPCWestPsiRaw%i",iHarm), Form("hTPCWestPsiRaw%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- // Rec
- hTPCFullQxRec[iHarm] = new TH2D(Form("hTPCFullQxRec%i",iHarm), Form("hTPCFullQxRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCFullQyRec[iHarm] = new TH2D(Form("hTPCFullQyRec%i",iHarm), Form("hTPCFullQyRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCEastQxRec[iHarm] = new TH2D(Form("hTPCEastQxRec%i",iHarm), Form("hTPCEastQxRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCEastQyRec[iHarm] = new TH2D(Form("hTPCEastQyRec%i",iHarm), Form("hTPCEastQyRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCWestQxRec[iHarm] = new TH2D(Form("hTPCWestQxRec%i",iHarm), Form("hTPCWestQxRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCWestQyRec[iHarm] = new TH2D(Form("hTPCWestQyRec%i",iHarm), Form("hTPCWestQyRec%i",iHarm), 200, -12,12, 16, -0.5, 15.5);
- hTPCFullPsiRec[iHarm] = new TH2D(Form("hTPCFullPsiRec%i",iHarm), Form("hTPCFullPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hTPCEastPsiRec[iHarm] = new TH2D(Form("hTPCEastPsiRec%i",iHarm), Form("hTPCEastPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hTPCWestPsiRec[iHarm] = new TH2D(Form("hTPCWestPsiRec%i",iHarm), Form("hTPCWestPsiRec%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- // Shift
- hTPCFullPsiShift[iHarm] = new TH2D(Form("hTPCFullPsiShift%i",iHarm), Form("hTPCFullPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hTPCEastPsiShift[iHarm] = new TH2D(Form("hTPCEastPsiShift%i",iHarm), Form("hTPCEastPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- hTPCWestPsiShift[iHarm] = new TH2D(Form("hTPCWestPsiShift%i",iHarm), Form("hTPCWestPsiShift%i",iHarm), 720, 0., 4.*TMath::Pi()/(iHarm+1.), 16, -0.5, 15.5);
- }
- for (int i = 0; i < 2; i++)
- {
- //Eventwise
- 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.);
- hVtxZ[i] = new TH1D(Form("hVtxZ%i", i), Form("vertex Z %i;Vtx_{Z}, [cm];N_{counts}", i), 200, -40., 40.);
- hVpdZ[i] = new TH1D(Form("hVpdZ%i", i), Form("vertex Z %i;Vpd_{Z}, [cm];N_{counts}", i), 200, -40., 40.);
- hVtxVpdZ[i] = new TH1D(Form("hVtxVpdZ%i", i), Form("VtxZ - VpdVz %i;Vtx_{Z}-vpdV_{Z}, [cm];N_{counts}", i), 200, -5., 5.);
- pVtxZvsRun[i] = new TProfile(Form("pVtxZvsRun%i",i),Form("pVtxZvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pVtxXvsRun[i] = new TProfile(Form("pVtxXvsRun%i",i),Form("pVtxXvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pVtxYvsRun[i] = new TProfile(Form("pVtxYvsRun%i",i),Form("pVtxYvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pVpdZvsRun[i] = new TProfile(Form("pVpdZvsRun%i",i),Form("pVpdZvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pVtxZVpdZvsRun[i] = new TProfile(Form("pVtxZVpdZvsRun%i",i),Form("pVtxZVpdZvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pRefMultvsRun[i] = new TProfile(Form("pRefMultvsRun%i",i),Form("pRefMultvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pAdcvsRun[i] = new TProfile(Form("pAdcvsRun%i",i),Form("pAdcvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pTofMultvsRun[i] = new TProfile(Form("pTofMultvsRun%i",i),Form("pTofMultvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pTofMatchvsRun[i] = new TProfile(Form("pTofMatchvsRun%i",i),Form("pTofMatchvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
-
- //Trackwise
- hNhitsFit[i] = new TH1D(Form("hNhitsFit%i", i), Form("N_{hist}^{Fit} %i;N_{hist}^{Fit};N_{counts}", i), 50, 0., 50.);
- hNhitsPoss[i] = new TH1D(Form("hNhitsPoss%i", i), Form("N_{hist}^{Poss} %i;N_{hist}^{Poss};N_{counts}", i), 50, 0., 50.);
- 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.);
- hPt[i] = new TH1D(Form("hPt%i", i), Form("p_{T} %i;p_{T}, [GeV/c];N_{counts}", i), 600, 0., 6.);
- hPtot[i] = new TH1D(Form("hPtot%i", i), Form("p_{tot} %i;p, [GeV/c];N_{counts}", i), 600, 0., 12.);
- hPhi[i] = new TH1D(Form("hPhi%i", i), Form("#varphi %i;#varphi, [rad];N_{counts}", i), 600, 0., 6.);
- hEta[i] = new TH1D(Form("hEta%i", i), Form("#eta %i;#eta;N_{counts}", i), 400, -2., 2.);
- hDCA[i] = new TH1D(Form("hDCA%i", i), Form("|DCA| %i;dca, [cm];N_{counts}", i), 200, 0., 5.);
- pPtvsRun[i] = new TProfile(Form("pPtvsRun%i",i),Form("pPtvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pEtavsRun[i] = new TProfile(Form("pEtavsRun%i",i),Form("pEtavsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pPhivsRun[i] = new TProfile(Form("pPhivsRun%i",i),Form("pPhivsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pNhitsFitvsRun[i] = new TProfile(Form("pNhitsFitvsRun%i",i),Form("pNhitsFitvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pNhitsPossvsRun[i] = new TProfile(Form("pNhitsPossvsRun%i",i),Form("pNhitsPossvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pDCAvsRun[i] = new TProfile(Form("pDCAvsRun%i",i),Form("pDCAvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pdEdxvsRun[i] = new TProfile(Form("pdEdxvsRun%i",i),Form("pdEdxvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pNtrPrimvsRun[i] = new TProfile(Form("pNtrPrimvsRun%i",i),Form("pNtrPrimvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pNtrGlobvsRun[i] = new TProfile(Form("pNtrGlobvsRun%i",i),Form("pNtrGlobvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pQx2TPCvsRun[i] = new TProfile(Form("pQx2TPCvsRun%i",i),Form("pQx2TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pQy2TPCvsRun[i] = new TProfile(Form("pQy2TPCvsRun%i",i),Form("pQy2TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pQx3TPCvsRun[i] = new TProfile(Form("pQx3TPCvsRun%i",i),Form("pQx3TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pQy3TPCvsRun[i] = new TProfile(Form("pQy3TPCvsRun%i",i),Form("pQy3TPCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pQx2BBCvsRun[i] = new TProfile(Form("pQx2BBCvsRun%i",i),Form("pQx2BBCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
- pQy2BBCvsRun[i] = new TProfile(Form("pQy2BBCvsRun%i",i),Form("pQy2BBCvsRun%i",i),fNBins, runId_min - 0.5, runId_max + 0.5);
-
- for (int iHarm=0; iHarm<2*fNharmonics; iHarm++)
- {
- 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);
- 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);
- 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);
- 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);
- 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);
- 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);
- }
- 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);
- 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);
- //PID
- 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);
- 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);
- for (int iPID = 0; iPID < 4; iPID++)
- {
- 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);
- 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);
- 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.);
- 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);
- 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.);
- }
- }
- for (int iPID = 0; iPID < 4; iPID++)
- {
- for (int iCent=0;iCent<16;iCent++)
- {
- 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);
- // hist[iPID][cent]->Fill(run,pt);
- }
- }
- //Counters
- Int_t fUsedTracks = 0;
- Int_t fQcountFull = 0;
- Int_t fQcountFullEast = 0;
- Int_t fQcountFullWest = 0;
- Int_t fQcountWest[NEtaGaps];
- Int_t fQcountEast[NEtaGaps];
- //Q-vectors
- TVector2 fQv2FullEP;
- TVector2 fQv2FullEPRaw;
- TVector2 fQv2FullSP;
- TVector2 fQv3FullEP;
- TVector2 fQv3FullEPRaw;
- TVector2 fQv3FullSP;
- TVector2 fQv2WestEP[NEtaGaps];
- TVector2 fQv2EastEP[NEtaGaps];
- TVector2 fQv2WestEPRaw[NEtaGaps];
- TVector2 fQv2EastEPRaw[NEtaGaps];
- TVector2 fQv2WestSP[NEtaGaps];
- TVector2 fQv2EastSP[NEtaGaps];
- TVector2 fQv3WestEP[NEtaGaps];
- TVector2 fQv3EastEP[NEtaGaps];
- TVector2 fQv3WestEPRaw[NEtaGaps];
- TVector2 fQv3EastEPRaw[NEtaGaps];
- TVector2 fQv3WestSP[NEtaGaps];
- TVector2 fQv3EastSP[NEtaGaps];
- TVector2 QvTMP;
- TVector2 fQv2WestSP_Shift[NEtaGaps];
- TVector2 fQv2EastSP_Shift[NEtaGaps];
- TVector2 fQv3WestSP_Shift[NEtaGaps];
- TVector2 fQv3EastSP_Shift[NEtaGaps];
- TVector2 fUvSP;
- Double_t weight;
- Int_t VtxSign;
- Double_t psiShiftedTMP;
- Double_t PsiSin;
- Double_t PsiCos;
- Long64_t goodEventCounter = 0;
- Double_t Psi2FullEP, Psi3FullEP;
- Double_t Psi2EastEP[NEtaGaps], Psi3EastEP[NEtaGaps];
- Double_t Psi2WestEP[NEtaGaps], Psi3WestEP[NEtaGaps];
- Double_t Psi2EastSP[NEtaGaps], Psi3EastSP[NEtaGaps];
- Double_t Psi2WestSP[NEtaGaps], Psi3WestSP[NEtaGaps];
- Double_t Psi_ReCenter;
- Int_t bin;
- Double_t res2EP, res3EP, res2EPsqr, res3EPsqr;
- Double_t flow2EP;
- Double_t flow3EP;
- Double_t res2SP, res3SP, res2SPsqr, res3SPsqr;
- Double_t flow2SP;
- Double_t flow3SP;
- Double_t psiSP_rec;
- int k = 0;
- TVector2 qvBBC_East, qvBBC_West;
- Int_t adcBin, shiftBin;
- Double_t psiEP_rec, psiEP_BBC_Full[fNharmonics+1], psiEP_BBC_East[fNharmonics+1], psiEP_BBC_West[fNharmonics+1], shiftCos, shiftSin;
- Double_t dPsi;
- // Loop over events
- for (Long64_t iEvent = 0; iEvent < events2read; iEvent++)
- {
- if ((iEvent + 1) % 1000 == 0)
- {
- std::cout << "Working on event #[" << (iEvent + 1)
- << "/" << events2read << "]" << std::endl;
- }
- Bool_t readEvent = femtoReader->readFemtoEvent(iEvent);
- if (!readEvent)
- {
- std::cout << "Something went wrong, Master! Nothing to analyze..." << std::endl;
- break;
- }
- // Retrieve femtoDst
- StFemtoDst *dst = femtoReader->femtoDst();
- // Retrieve event information
- StFemtoEvent *event = dst->event();
- if (!event)
- {
- std::cout << "Something went wrong, Master! Event is hiding from me..." << std::endl;
- break;
- }
- TVector3 pVtx = event->primaryVertex();
- hVtxXY[0]->Fill(pVtx.X(), pVtx.Y());
- hVtxZ[0]->Fill(pVtx.Z());
- hVtxVpdZ[0]->Fill(pVtx.Z() - event->vpdVz());
- hVpdZ[0]->Fill(event->vpdVz());
- pVtxXvsRun[0]->Fill(event->runId(),pVtx.X());
- pVtxYvsRun[0]->Fill(event->runId(),pVtx.Y());
- pVtxZvsRun[0]->Fill(event->runId(),pVtx.Z());
- pVpdZvsRun[0]->Fill(event->runId(),event->vpdVz());
- pVtxZVpdZvsRun[0]->Fill(event->runId(),pVtx.Z() - event->vpdVz());
- pRefMultvsRun[0]->Fill(event->runId(),event->refMult());
- pTofMultvsRun[0]->Fill(event->runId(),event->numberOfBTofHit());
- pTofMatchvsRun[0]->Fill(event->runId(),event->numberOfTofMatched());
- pNtrPrimvsRun[0]->Fill(event->runId(),event->numberOfPrimaryTracks());
- pNtrGlobvsRun[0]->Fill(event->runId(),event->numberOfGlobalTracks());
- Int_t fAdcSum=0;
- for (int i=0;i<16;i++)
- {
- fAdcSum += event->bbcAdcEast(i);
- fAdcSum += event->bbcAdcWest(i);
- }
- pAdcvsRun[0]->Fill(event->runId(),fAdcSum);
- // Event selection
- if (!isGoodEvent(event))
- continue;
- VtxSign = GetVzBin(event->primaryVertex().Z());
- if (VtxSign == -1)
- continue;
- hNevVsRun->Fill(event->runId(),event->cent16());
- //Counters
- fUsedTracks = 0;
- fQcountFull = 0;
- fQcountFullEast = 0;
- fQcountFullWest = 0;
- psiShiftedTMP = 0.;
- PsiSin = 0.;
- PsiCos = 0.;
- fQv2FullEP.Set(0.,0.);
- fQv2FullEPRaw.Set(0.,0.);
- fQv2FullSP.Set(0.,0.);
- fQv3FullEP.Set(0.,0.);
- fQv3FullEPRaw.Set(0.,0.);
- fQv3FullSP.Set(0.,0.);
- QvTMP.Set(0.,0.);
- for (int i=0; i<NEtaGaps; i++)
- {
- fQcountWest[i] = 0;
- fQcountEast[i] = 0;
- fQv2WestEP[i].Set(0.,0.);
- fQv2EastEP[i].Set(0.,0.);
- fQv2WestEPRaw[i].Set(0.,0.);
- fQv2EastEPRaw[i].Set(0.,0.);
- fQv2WestSP[i].Set(0.,0.);
- fQv2EastSP[i].Set(0.,0.);
- fQv3WestEP[i].Set(0.,0.);
- fQv3EastEP[i].Set(0.,0.);
- fQv3WestEPRaw[i].Set(0.,0.);
- fQv3EastEPRaw[i].Set(0.,0.);
- fQv3WestSP[i].Set(0.,0.);
- fQv3EastSP[i].Set(0.,0.);
- }
- hVtxXY[1]->Fill(pVtx.X(), pVtx.Y());
- hVtxZ[1]->Fill(pVtx.Z());
- hVtxVpdZ[1]->Fill(pVtx.Z() - event->vpdVz());
- hVpdZ[1]->Fill(event->vpdVz());
- pVtxXvsRun[1]->Fill(event->runId(),pVtx.X());
- pVtxYvsRun[1]->Fill(event->runId(),pVtx.Y());
- pVtxZvsRun[1]->Fill(event->runId(),pVtx.Z());
- pVpdZvsRun[1]->Fill(event->runId(),event->vpdVz());
- pVtxZVpdZvsRun[1]->Fill(event->runId(),pVtx.Z() - event->vpdVz());
- pRefMultvsRun[1]->Fill(event->runId(),event->refMult());
- pTofMultvsRun[1]->Fill(event->runId(),event->numberOfBTofHit());
- pTofMatchvsRun[1]->Fill(event->runId(),event->numberOfTofMatched());
- pNtrPrimvsRun[1]->Fill(event->runId(),event->numberOfPrimaryTracks());
- pNtrGlobvsRun[1]->Fill(event->runId(),event->numberOfGlobalTracks());
- fAdcSum=0;
- for (int i=0;i<16;i++)
- {
- fAdcSum += event->bbcAdcEast(i);
- fAdcSum += event->bbcAdcWest(i);
- }
- pAdcvsRun[1]->Fill(event->runId(),fAdcSum);
- if (k == 0)
- {
- std::cout << "East BBC tileId = 8 with phi = " << BBC_GetPhi(0, 8) << std::endl;
- std::cout << "West BBC tileId = 8 with phi = " << BBC_GetPhi(1, 8) << std::endl;
- k++;
- }
- //BBC EP
- //Raw
- TVector2 qRawEast[fNharmonics+1], qRawWest[fNharmonics+1], qRawFull[fNharmonics+1];
- Double_t psiRawEastBBC[fNharmonics+1], psiRawWestBBC[fNharmonics+1], psiRawFullBBC[fNharmonics+1];
- float nHitE_Raw[16], nHitW_Raw[16];
- // 0 is for full event plane, 1 is for west only event plane, 2 is for east only event plane
- //the number of hits in the bbc tile is simply the adc
- for(int i=0;i<16;i++){
- nHitE_Raw[i] = event->bbcAdcEast(i);
- nHitW_Raw[i] = event->bbcAdcWest(i);
- }
- for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
- {
- //BBC Full
- //v1*y > 0 by conventions most reasonable to the STAR event plane
- //therefore the west event plane has the right sign for the event plane
- //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
- //The east event plane points opposite that of the west and the west finds the "correct" EP for the STAR coordinates, use West-East
- Double_t QVector_RawX = 0., QVector_RawY = 0.;
- Float_t eXFsum_Raw = 0., wXFsum_Raw = 0., eYFsum_Raw = 0., wYFsum_Raw = 0., eWFgt_Raw=0., wWFgt_Raw = 0.;
- TVector2 QFullRaw;
- for(int iTile = 0; iTile < 16; iTile++) {
- eXFsum_Raw += cos((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
- wXFsum_Raw += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
- wYFsum_Raw += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
- eYFsum_Raw += sin((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
- eWFgt_Raw += nHitE_Raw[iTile];
- wWFgt_Raw += nHitW_Raw[iTile];
- }
- QVector_RawX=(eWFgt_Raw > 0. && wWFgt_Raw > 0.) ? wXFsum_Raw/wWFgt_Raw - eXFsum_Raw/eWFgt_Raw:0.;
- QVector_RawY=(eWFgt_Raw > 0. && wWFgt_Raw > 0.) ? wYFsum_Raw/wWFgt_Raw - eYFsum_Raw/eWFgt_Raw:0.;
- QFullRaw.Set(QVector_RawX,QVector_RawY);
- qRawFull[iHarm] = QFullRaw;
- if(qRawFull[iHarm].Mod()){
- psiRawFullBBC[iHarm] = qRawFull[iHarm].Phi()/(iHarm+1.);
- if(psiRawFullBBC[iHarm] < 0.0) {psiRawFullBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- //BBC West
- //v1*y > 0 by conventions most reasonable to the STAR event plane
- //therefore the west event plane has the right sign for the event plane
- //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
- Double_t QVectorW_RawX = 0., QVectorW_RawY = 0.;
- Float_t wXsum_Raw = 0., wYsum_Raw = 0., wWgt_Raw = 0.;
- Double_t EventPlaneRawWest = 0.;
- TVector2 QWestRaw;
- for(int iTile = 0; iTile < 16; iTile++) {
- wXsum_Raw += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
- wYsum_Raw += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Raw[iTile];
- wWgt_Raw += nHitW_Raw[iTile];
- }
- QVectorW_RawX = (wWgt_Raw > 0.0) ? wXsum_Raw/wWgt_Raw:0.0;
- QVectorW_RawY = (wWgt_Raw > 0.0) ? wYsum_Raw/wWgt_Raw:0.0;
-
- QWestRaw.Set(QVectorW_RawX,QVectorW_RawY);
- qRawWest[iHarm] = QWestRaw;
- if(qRawWest[iHarm].Mod()){
- psiRawWestBBC[iHarm] = qRawWest[iHarm].Phi()/(iHarm+1.);
- if(psiRawWestBBC[iHarm] < 0.0) {psiRawWestBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- //BBC East
- //in BBC_GetPhi the fact that the tile numbering scheme is reversed about the y axis for the east BBC is accounted for
- //The east event plane points opposite that of the west since the spectators are on different sides -- see west comments
- Double_t QVectorE_RawX = 0., QVectorE_RawY = 0.;
- Float_t eXsum_Raw = 0., eYsum_Raw = 0., eWgt_Raw = 0.;
- TVector2 QEastRaw;
- for(int iTile = 0; iTile < 16; iTile++) {
- eXsum_Raw += cos(BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
- eYsum_Raw += sin(BBC_GetPhi(0,iTile))*nHitE_Raw[iTile];
- eWgt_Raw += nHitE_Raw[iTile];
- }
- QVectorE_RawX = (eWgt_Raw > 0.0) ? eXsum_Raw/eWgt_Raw:0.0;
- QVectorE_RawY = (eWgt_Raw > 0.0) ? eYsum_Raw/eWgt_Raw:0.0;
-
- QEastRaw.Set(QVectorE_RawX,QVectorE_RawY);
- qRawEast[iHarm] = QEastRaw;
- if(qRawEast[iHarm].Mod()){
- psiRawEastBBC[iHarm] = qRawEast[iHarm].Phi()/(iHarm+1.);
- if(psiRawEastBBC[iHarm] < 0.0) {psiRawEastBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- }
- for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
- {
- hBBCFullQxRaw[iHarm]->Fill(qRawFull[iHarm].X(),event->cent16());
- hBBCFullQyRaw[iHarm]->Fill(qRawFull[iHarm].Y(),event->cent16());
- hBBCEastQxRaw[iHarm]->Fill(qRawEast[iHarm].X(),event->cent16());
- hBBCEastQyRaw[iHarm]->Fill(qRawEast[iHarm].Y(),event->cent16());
- hBBCWestQxRaw[iHarm]->Fill(qRawWest[iHarm].X(),event->cent16());
- hBBCWestQyRaw[iHarm]->Fill(qRawWest[iHarm].Y(),event->cent16());
- hBBCFullPsiRaw[iHarm]->Fill(psiRawFullBBC[iHarm], event->cent16());
- hBBCEastPsiRaw[iHarm]->Fill(psiRawEastBBC[iHarm], event->cent16());
- hBBCWestPsiRaw[iHarm]->Fill(psiRawWestBBC[iHarm], event->cent16());
- }
- // Gain & Rec & Shift
- Double_t psiGainEastBBC[fNharmonics+1], psiGainWestBBC[fNharmonics+1], psiGainFullBBC[fNharmonics+1];
- Double_t psiRecEastBBC[fNharmonics+1], psiRecWestBBC[fNharmonics+1], psiRecFullBBC[fNharmonics+1];
- Double_t adcNormEastInner = BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),1,6)/6.;
- Double_t adcNormEastOuter = BBCgain_East[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),7,16)/10.;
- Double_t adcNormWestInner = BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),1,6)/6.;
- Double_t adcNormWestOuter = BBCgain_West[VtxSign]->Integral(BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),BBCgain_East[VtxSign]->GetXaxis()->FindBin(event->runId()),7,16)/10.;
- float egain[16]={
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),1))/adcNormEastInner),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),2))/adcNormEastInner),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),3))/adcNormEastInner),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),4))/adcNormEastInner),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),5))/adcNormEastInner),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),6))/adcNormEastInner),
- (float)((0.5)*BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),7))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),8))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),9))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),10))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),11))/adcNormEastOuter),
- (float)((0.5)*BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),12))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),13))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),14))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),15))/adcNormEastOuter),
- (float)(BBCgain_East[VtxSign]->GetBinContent(BBCgain_East[VtxSign]->FindBin(event->runId(),16))/adcNormEastOuter)
- };
- float wgain[16]={
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),1))/adcNormWestInner),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),2))/adcNormWestInner),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),3))/adcNormWestInner),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),4))/adcNormWestInner),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),5))/adcNormWestInner),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),6))/adcNormWestInner),
- (float)((0.5)*BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),7))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),8))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),9))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),10))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),11))/adcNormWestOuter),
- (float)((0.5)*BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),12))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),13))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),14))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),15))/adcNormWestOuter),
- (float)(BBCgain_West[VtxSign]->GetBinContent(BBCgain_West[VtxSign]->FindBin(event->runId(),16))/adcNormWestOuter)
- };
- float nHitE_Gain[16], nHitW_Gain[16];
- for(int i=0;i<16;i++){
- nHitE_Gain[i] = event->bbcAdcEast(i)/egain[i];
- nHitW_Gain[i] = event->bbcAdcWest(i)/wgain[i];
- }
- Double_t psiWest[fNharmonics+1];
- Double_t psiEast[fNharmonics+1];
- Double_t psiFull[fNharmonics+1];
- TVector2 qGainEast[fNharmonics+1], qGainWest[fNharmonics+1], qGainFull[fNharmonics+1];
- TVector2 qRecEast[fNharmonics+1], qRecWest[fNharmonics+1], qRecFull[fNharmonics+1];
- for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
- {
- //BBC Full
- //v1*y > 0 by conventions most reasonable to the STAR event plane
- //therefore the west event plane has the right sign for the event plane
- //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
- //The east event plane points opposite that of the west and the west finds the "correct" EP for the STAR coordinates, use West-East
- Double_t QVector_GainX = 0., QVector_GainY = 0.;
- Double_t EventPlaneGainFull = 0.;
- Float_t eXFsum_Gain = 0., wXFsum_Gain = 0., eYFsum_Gain = 0., wYFsum_Gain = 0., eWFgt_Gain=0., wWFgt_Gain = 0.;
- TVector2 QFullGain;
- for(int iTile = 0; iTile < 16; iTile++) {
- eXFsum_Gain += cos((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
- wXFsum_Gain += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
- eYFsum_Gain += sin((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
- wYFsum_Gain += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
- eWFgt_Gain += nHitE_Gain[iTile];
- wWFgt_Gain += nHitW_Gain[iTile];
- }
- QVector_GainX=(eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wXFsum_Gain/wWFgt_Gain + eXFsum_Gain/eWFgt_Gain:0.;
- QVector_GainY=(eWFgt_Gain > 0. && wWFgt_Gain > 0.) ? wYFsum_Gain/wWFgt_Gain + eYFsum_Gain/eWFgt_Gain:0.;
- QFullGain.Set(QVector_GainX,QVector_GainY);
- if(QFullGain.Mod()){
- EventPlaneGainFull = QFullGain.Phi();
- if(EventPlaneGainFull < 0.0) {EventPlaneGainFull += 2.*TMath::Pi()/(iHarm+1.);}
- }
- psiFull[iHarm] = EventPlaneGainFull;
- qGainFull[iHarm] = QFullGain;
- if(qGainFull[iHarm].Mod()){
- psiGainFullBBC[iHarm] = qGainFull[iHarm].Phi()/(iHarm+1.);
- if(psiGainFullBBC[iHarm] < 0.0) {psiGainFullBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- //BBC West
- //v1*y > 0 by conventions most reasonable to the STAR event plane
- //therefore the west event plane has the right sign for the event plane
- //more simply (Psi_w = 0 => QWest . x_STAR > 0) and (Psi_e = 0 => QEast . x_STAR < 0)
- Double_t QVectorW_GainX = 0., QVectorW_GainY = 0.;
- Double_t EventPlaneGainWest = 0.;
- Float_t wXsum_Gain = 0., wYsum_Gain = 0., wWgt_Gain = 0.;
- TVector2 QWestGain;
- for(int iTile = 0; iTile < 16; iTile++) {
- wXsum_Gain += cos((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
- wYsum_Gain += sin((iHarm+1.)*BBC_GetPhi(1,iTile))*nHitW_Gain[iTile];
- wWgt_Gain += nHitW_Gain[iTile];
- }
- QVectorW_GainX = (wWgt_Gain > 0.0) ? wXsum_Gain/wWgt_Gain:0.0;
- QVectorW_GainY = (wWgt_Gain > 0.0) ? wYsum_Gain/wWgt_Gain:0.0;
-
- QWestGain.Set(QVectorW_GainX,QVectorW_GainY);
- if(QWestGain.Mod()){
- EventPlaneGainWest = QWestGain.Phi()/(iHarm+1.);
- if(EventPlaneGainWest < 0.0) {EventPlaneGainWest += 2.*TMath::Pi()/(iHarm+1.);}
- }
- psiWest[iHarm] = EventPlaneGainWest;
- qGainWest[iHarm] = QWestGain;
- if(qGainWest[iHarm].Mod()){
- psiGainWestBBC[iHarm] = qGainWest[iHarm].Phi()/(iHarm+1.);
- if(psiGainWestBBC[iHarm] < 0.0) {psiGainWestBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- //BBC East
- //in BBC_GetPhi the fact that the tile numbering scheme is reversed about the y axis for the east BBC is accounted for
- //The east event plane points opposite that of the west since the spectators are on different sides -- see west comments
- Double_t QVectorE_GainX = 0., QVectorE_GainY = 0.;
- Double_t EventPlaneGainEast = 0.;
- Float_t eXsum_Gain = 0., eYsum_Gain = 0., eWgt_Gain = 0.;
- TVector2 QEastGain;
- for(int iTile = 0; iTile < 16; iTile++) {
- eXsum_Gain += cos((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
- eYsum_Gain += sin((iHarm+1.)*BBC_GetPhi(0,iTile))*nHitE_Gain[iTile];
- eWgt_Gain += nHitE_Gain[iTile];
- }
- QVectorE_GainX = (eWgt_Gain > 0.0) ? eXsum_Gain/eWgt_Gain:0.0;
- QVectorE_GainY = (eWgt_Gain > 0.0) ? eYsum_Gain/eWgt_Gain:0.0;
-
- QEastGain.Set(QVectorE_GainX,QVectorE_GainY);
- if(QEastGain.Mod()){
- EventPlaneGainEast = QEastGain.Phi()/(iHarm+1.);
- if(EventPlaneGainEast < 0.0) {EventPlaneGainEast += 2.*TMath::Pi()/(iHarm+1.);}
- }
- psiEast[iHarm] = EventPlaneGainEast;
- qGainEast[iHarm] = QEastGain;
- if(qGainEast[iHarm].Mod()){
- psiGainEastBBC[iHarm] = qGainEast[iHarm].Phi()/(iHarm+1.);
- if(psiGainEastBBC[iHarm] < 0.0) {psiGainEastBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- }
- // Apply recentering
- for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
- {
- TVector2 recFullTMP((Double_t)p_BBC_Qx_Full_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
- (Double_t)p_BBC_Qy_Full_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- TVector2 recEastTMP((Double_t)p_BBC_Qx_East_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
- (Double_t)p_BBC_Qy_East_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- TVector2 recWestTMP((Double_t)p_BBC_Qx_West_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()],
- (Double_t)p_BBC_Qy_West_EP[iHarm][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- qRecFull[iHarm] = qGainFull[iHarm] - recFullTMP;
- qRecEast[iHarm] = qGainEast[iHarm] - recEastTMP;
- qRecWest[iHarm] = qGainWest[iHarm] - recWestTMP;
- if(qRecFull[iHarm].Mod()){
- psiRecFullBBC[iHarm] = qRecFull[iHarm].Phi()/(iHarm+1.);
- if(psiRecFullBBC[iHarm] < 0.0) {psiRecFullBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- if(qRecEast[iHarm].Mod()){
- psiRecEastBBC[iHarm] = qRecEast[iHarm].Phi()/(iHarm+1.);
- if(psiRecEastBBC[iHarm] < 0.0) {psiRecEastBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- if(qRecWest[iHarm].Mod()){
- psiRecWestBBC[iHarm] = qRecWest[iHarm].Phi()/(iHarm+1.);
- if(psiRecWestBBC[iHarm] < 0.0) {psiRecWestBBC[iHarm] += 2.*TMath::Pi()/(iHarm+1.);}
- }
- }
- for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- {
- // Full BBC
- if (qRecFull[iHarm].Mod())
- {
- psiEP_rec = qRecFull[iHarm].Phi() / (iHarm+1.);
- if (psiEP_rec < 0.) {psiEP_rec += 2*TMath::Pi()/(iHarm+1.);}
- }
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax*3; iK++)
- {
- shiftCos = p_BBC_Cos_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- shiftSin = p_BBC_Sin_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- dPsi += (1. / (iHarm+1.)) * (2. / (iK + 1.)) *
- (-1. * shiftSin * cos((iHarm+1.) * (iK + 1.) * psiEP_rec) +
- 1. * shiftCos * sin((iHarm+1.) * (iK + 1.) * psiEP_rec));
- }
- psiEP_BBC_Full[iHarm] = psiEP_rec + dPsi;
- // East BBC
- if (qRecEast[iHarm].Mod())
- {
- psiEP_rec = qRecEast[iHarm].Phi() / (iHarm+1.);
- if (psiEP_rec < 0.) {psiEP_rec += 2*TMath::Pi()/(iHarm+1.);}
- }
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax*3; iK++)
- {
- shiftCos = p_BBC_Cos_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- shiftSin = p_BBC_Sin_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- dPsi += (1. / (iHarm+1.)) * (2. / (iK + 1.)) *
- (-1. * shiftSin * cos((iHarm+1.) * (iK + 1.) * psiEP_rec) +
- 1. * shiftCos * sin((iHarm+1.) * (iK + 1.) * psiEP_rec));
- }
- psiEP_BBC_East[iHarm] = psiEP_rec + dPsi;
- // West BBC
- if (qRecWest[iHarm].Mod())
- {
- psiEP_rec = qRecWest[iHarm].Phi() / (iHarm+1.);
- if (psiEP_rec < 0.) {psiEP_rec += 2*TMath::Pi()/(iHarm+1.);}
- }
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax*3; iK++)
- {
- shiftCos = p_BBC_Cos_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- shiftSin = p_BBC_Sin_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- dPsi += (1. / (iHarm+1.)) * (2. / (iK + 1.)) *
- (-1. * shiftSin * cos((iHarm+1.) * (iK + 1.) * psiEP_rec) +
- 1. * shiftCos * sin((iHarm+1.) * (iK + 1.) * psiEP_rec));
- }
- psiEP_BBC_West[iHarm] = psiEP_rec + dPsi;
- if(psiEP_BBC_Full[iHarm] < 0.0) psiEP_BBC_Full[iHarm] += 2.*TMath::Pi()/(iHarm+1.);
- if(psiEP_BBC_Full[iHarm] > 2.*TMath::Pi()/(iHarm+1.)) psiEP_BBC_Full[iHarm] -= 2.*TMath::Pi()/(iHarm+1.);
- if(psiEP_BBC_East[iHarm] < 0.0) psiEP_BBC_East[iHarm] += 2.*TMath::Pi()/(iHarm+1.);
- if(psiEP_BBC_East[iHarm] > 2.*TMath::Pi()/(iHarm+1.)) psiEP_BBC_East[iHarm] -= 2.*TMath::Pi()/(iHarm+1.);
- if(psiEP_BBC_West[iHarm] < 0.0) psiEP_BBC_West[iHarm] += 2.*TMath::Pi()/(iHarm+1.);
- if(psiEP_BBC_West[iHarm] > 2.*TMath::Pi()/(iHarm+1.)) psiEP_BBC_West[iHarm] -= 2.*TMath::Pi()/(iHarm+1.);
- }
- for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
- {
- hBBCFullQxGain[iHarm]->Fill(qGainFull[iHarm].X(),event->cent16());
- hBBCFullQyGain[iHarm]->Fill(qGainFull[iHarm].Y(),event->cent16());
- hBBCEastQxGain[iHarm]->Fill(qGainEast[iHarm].X(),event->cent16());
- hBBCEastQyGain[iHarm]->Fill(qGainEast[iHarm].Y(),event->cent16());
- hBBCWestQxGain[iHarm]->Fill(qGainWest[iHarm].X(),event->cent16());
- hBBCWestQyGain[iHarm]->Fill(qGainWest[iHarm].Y(),event->cent16());
- hBBCFullPsiGain[iHarm]->Fill(psiGainFullBBC[iHarm], event->cent16());
- hBBCEastPsiGain[iHarm]->Fill(psiGainEastBBC[iHarm], event->cent16());
- hBBCWestPsiGain[iHarm]->Fill(psiGainWestBBC[iHarm], event->cent16());
- hBBCFullQxRec[iHarm]->Fill(qRecFull[iHarm].X(),event->cent16());
- hBBCFullQyRec[iHarm]->Fill(qRecFull[iHarm].Y(),event->cent16());
- hBBCEastQxRec[iHarm]->Fill(qRecEast[iHarm].X(),event->cent16());
- hBBCEastQyRec[iHarm]->Fill(qRecEast[iHarm].Y(),event->cent16());
- hBBCWestQxRec[iHarm]->Fill(qRecWest[iHarm].X(),event->cent16());
- hBBCWestQyRec[iHarm]->Fill(qRecWest[iHarm].Y(),event->cent16());
- hBBCFullPsiRec[iHarm]->Fill(psiRecFullBBC[iHarm], event->cent16());
- hBBCEastPsiRec[iHarm]->Fill(psiRecEastBBC[iHarm], event->cent16());
- hBBCWestPsiRec[iHarm]->Fill(psiRecWestBBC[iHarm], event->cent16());
- hBBCFullPsiShift[iHarm]->Fill(psiEP_BBC_Full[iHarm], event->cent16());
- hBBCEastPsiShift[iHarm]->Fill(psiEP_BBC_East[iHarm], event->cent16());
- hBBCWestPsiShift[iHarm]->Fill(psiEP_BBC_West[iHarm], event->cent16());
-
- // East
- pQx2BBCvsRun[0]->Fill(event->runId(),qRecEast[1].X());
- pQy2BBCvsRun[0]->Fill(event->runId(),qRecEast[1].Y());
- pCosBBCvsRun[0][0]->Fill(event->runId(),TMath::Cos(1*psiEP_BBC_East[1]));
- pCosBBCvsRun[1][0]->Fill(event->runId(),TMath::Cos(2*psiEP_BBC_East[1]));
- pCosBBCvsRun[2][0]->Fill(event->runId(),TMath::Cos(3*psiEP_BBC_East[1]));
- pCosBBCvsRun[3][0]->Fill(event->runId(),TMath::Cos(4*psiEP_BBC_East[1]));
- pSinBBCvsRun[0][0]->Fill(event->runId(),TMath::Sin(1*psiEP_BBC_East[1]));
- pSinBBCvsRun[1][0]->Fill(event->runId(),TMath::Sin(2*psiEP_BBC_East[1]));
- pSinBBCvsRun[2][0]->Fill(event->runId(),TMath::Sin(3*psiEP_BBC_East[1]));
- pSinBBCvsRun[3][0]->Fill(event->runId(),TMath::Sin(4*psiEP_BBC_East[1]));
- // West
- pQx2BBCvsRun[1]->Fill(event->runId(),qRecWest[1].X());
- pQy2BBCvsRun[1]->Fill(event->runId(),qRecWest[1].Y());
- pCosBBCvsRun[0][1]->Fill(event->runId(),TMath::Cos(1*psiEP_BBC_West[1]));
- pCosBBCvsRun[1][1]->Fill(event->runId(),TMath::Cos(2*psiEP_BBC_West[1]));
- pCosBBCvsRun[2][1]->Fill(event->runId(),TMath::Cos(3*psiEP_BBC_West[1]));
- pCosBBCvsRun[3][1]->Fill(event->runId(),TMath::Cos(4*psiEP_BBC_West[1]));
- pSinBBCvsRun[0][1]->Fill(event->runId(),TMath::Sin(1*psiEP_BBC_West[1]));
- pSinBBCvsRun[1][1]->Fill(event->runId(),TMath::Sin(2*psiEP_BBC_West[1]));
- pSinBBCvsRun[2][1]->Fill(event->runId(),TMath::Sin(3*psiEP_BBC_West[1]));
- pSinBBCvsRun[3][1]->Fill(event->runId(),TMath::Sin(4*psiEP_BBC_West[1]));
- }
- // Track analysis
- Int_t nTracks = dst->numberOfTracks();
- // Track loop
- for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
- {
- // Retrieve i-th femto track
- StFemtoTrack *femtoTrack = dst->track(iTrk);
- if (!femtoTrack)
- continue;
- hNhitsFit[0]->Fill(femtoTrack->nHitsFit());
- hNhitsPoss[0]->Fill(femtoTrack->nHitsPoss());
- if (femtoTrack->nHitsPoss() != 0)
- {
- hNhitsRatio[0]->Fill((Double_t)(femtoTrack->nHitsFit() / femtoTrack->nHitsPoss()));
- }
- hPt[0]->Fill(femtoTrack->pt());
- hPtot[0]->Fill(femtoTrack->ptot());
- hEta[0]->Fill(femtoTrack->eta());
- hPhi[0]->Fill(femtoTrack->phi());
- hDCA[0]->Fill(femtoTrack->gDCA(pVtx).Mag());
- pPtvsRun[0]->Fill(event->runId(),femtoTrack->pt());
- pEtavsRun[0]->Fill(event->runId(),femtoTrack->eta());
- pPhivsRun[0]->Fill(event->runId(),femtoTrack->phi());
- pNhitsFitvsRun[0]->Fill(event->runId(),femtoTrack->nHitsFit());
- pNhitsPossvsRun[0]->Fill(event->runId(),femtoTrack->nHitsPoss());
- pDCAvsRun[0]->Fill(event->runId(),femtoTrack->gDCA(pVtx).Mag());
- pdEdxvsRun[0]->Fill(event->runId(),femtoTrack->dEdx());
- // Must be a primary track
- if (!femtoTrack->isPrimary())
- continue;
-
- //Track selection
- if (isGoodTrack(femtoTrack, energy, pVtx))
- {
- fUsedTracks++;
-
- //Fill recentered Q-vectors for all TPC (FULL)
- weight = GetWeight(femtoTrack);
-
- QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qx_FULL_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- fQv2FullEP += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_FULL_EP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_FULL_EP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- fQv3FullEP += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_FULL_SP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- fQv2FullSP += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_FULL_SP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_FULL_SP[1][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()]);
- fQv3FullSP += 1.*(CalcQvector(femtoTrack,3) - QvTMP);
- QvTMP.Set(0.,0.);
-
- fQcountFull++;
- if (femtoTrack->eta() >= 0.)
- {
- fQcountFullWest++;
- }
- else
- {
- fQcountFullEast++;
- }
- for (int iGap=0; iGap<NEtaGaps; iGap++)
- {
- //fill Qvectors from EAST arm
- if (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < 0.)//cutEtaGap[iGap])
- {
- QvTMP.SetX(p_Qx_EAST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_EAST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv2EastEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_EAST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_EAST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv3EastEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_EAST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv2EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_EAST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_EAST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv3EastSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
- QvTMP.Set(0.,0.);
- fQcountEast[iGap]++;
- }
- //fill Qvectors from WEST arm
- if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 0.)//cutEtaGap[iGap])
- {
- QvTMP.SetX(p_Qx_WEST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_WEST_EP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv2WestEP[iGap] += weight * (CalcQvector(femtoTrack, 2) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_WEST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_WEST_EP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv3WestEP[iGap] += weight * (CalcQvector(femtoTrack, 3) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_WEST_SP[0][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv2WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 2) - QvTMP);
- QvTMP.Set(0., 0.);
- QvTMP.SetX(p_Qx_WEST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- QvTMP.SetY(p_Qy_WEST_SP[1][VtxSign][(Int_t)event->cent16()][iGap][(Double_t)event->runId()]);
- fQv3WestSP[iGap] += 1. * (CalcQvector(femtoTrack, 3) - QvTMP);
- QvTMP.Set(0.,0.);
- fQcountWest[iGap]++;
- }
- }
- }
- //Track selection
- if (!isGoodTrackFlow(femtoTrack, energy, pVtx))
- continue;
- hNhitsFit[1]->Fill(femtoTrack->nHitsFit());
- hNhitsPoss[1]->Fill(femtoTrack->nHitsPoss());
- if (femtoTrack->nHitsPoss() != 0)
- {
- hNhitsRatio[1]->Fill((Double_t)(femtoTrack->nHitsFit() / femtoTrack->nHitsPoss()));
- }
- hPt[1]->Fill(femtoTrack->pt());
- hPtot[1]->Fill(femtoTrack->ptot());
- hPhi[1]->Fill(femtoTrack->phi());
- hEta[1]->Fill(femtoTrack->eta());
- hDCA[1]->Fill(femtoTrack->gDCA(pVtx).Mag());
- pPtvsRun[1]->Fill(event->runId(),femtoTrack->pt());
- pEtavsRun[1]->Fill(event->runId(),femtoTrack->eta());
- pPhivsRun[1]->Fill(event->runId(),femtoTrack->phi());
- pNhitsFitvsRun[1]->Fill(event->runId(),femtoTrack->nHitsFit());
- pNhitsPossvsRun[1]->Fill(event->runId(),femtoTrack->nHitsPoss());
- pDCAvsRun[1]->Fill(event->runId(),femtoTrack->gDCA(pVtx).Mag());
- pdEdxvsRun[1]->Fill(event->runId(),femtoTrack->dEdx());
- Int_t i_part = -1;
- //TPC-only
- if (!femtoTrack->isTofTrack())
- {
- //pion id
- if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.6 &&
- TMath::Abs(femtoTrack->nSigmaPion()) < 2)
- {
- i_part = 0;
- }
- // kaon id
- if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 0.5 &&
- TMath::Abs(femtoTrack->nSigmaKaon()) < 2)
- {
- i_part = 1;
- }
- // proton id
- if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 0.9 &&
- TMath::Abs(femtoTrack->nSigmaProton()) < 2)
- {
- i_part = 2;
- }
- //pion id
- if (femtoTrack->ptot() >= 0.6 && femtoTrack->ptot() < 0.7 &&
- TMath::Abs(femtoTrack->nSigmaPion()) < 2 &&
- TMath::Abs(femtoTrack->nSigmaKaon()) > 2)
- {
- i_part = 0;
- }
- // kaon id
- if (femtoTrack->ptot() >= 0.5 && femtoTrack->ptot() < 0.7 &&
- TMath::Abs(femtoTrack->nSigmaKaon()) < 2 &&
- TMath::Abs(femtoTrack->nSigmaPion()) > 3)
- {
- i_part = 1;
- }
- // proton id
- if (femtoTrack->ptot() >= 0.9 && femtoTrack->ptot() < 1.2 &&
- TMath::Abs(femtoTrack->nSigmaProton()) < 2 &&
- TMath::Abs(femtoTrack->nSigmaPion()) > 3)
- {
- i_part = 2;
- }
- }
- //TPC+TOF
- if (isGoodPID(femtoTrack))
- {
- // pion id
- if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
- TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
- femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)
- {
- i_part = 0;
- }
- // kaon id
- if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
- TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
- femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)
- {
- i_part = 1;
- }
- // proton id
- if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
- TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
- femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)
- {
- i_part = 2;
- }
- }
- for (int iGap=0; iGap<NEtaGaps; iGap++)
- {
- //EAST OR WEST
- if ( (femtoTrack->eta() > -1.*cutEta && femtoTrack->eta() < -1.*cutEtaGap.at(iGap)) ||
- (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 1.*cutEtaGap.at(iGap)) )
- {
- //Fill recentered Q-vectors for all TPC (FULL)
- weight = GetWeight(femtoTrack);
- h_pt_spectra[3][(Int_t)event->cent16()]->Fill((Double_t)event->runId(),femtoTrack->pt(),weight);
- if (i_part != -1)
- {
- h_pt_spectra[i_part][(Int_t)event->cent16()]->Fill((Double_t)event->runId(),femtoTrack->pt(),weight);
- }
- }
- }
- // Check if track has TOF signal
- if (femtoTrack->isTofTrack())
- {
- if (femtoTrack->gDCA(pVtx).Mag() >= 1.)
- continue;
- if (femtoTrack->pt() > 0.5 && femtoTrack->pt() < 0.7)
- {
- hnSigM2_lowpt[0]->Fill(femtoTrack->nSigmaPion(), femtoTrack->massSqr());
- }
- if (femtoTrack->pt() > 2.2 && femtoTrack->pt() < 2.4)
- {
- hnSigM2_highpt[0]->Fill(femtoTrack->nSigmaPion(), femtoTrack->massSqr());
- }
- hdEdxQp[0][0]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
- hM2Qp[0][0]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
- hNsigQp[0][0]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
- hM2pt[0][0]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
- hNsigpt[0][0]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
- // pion
- if (femtoTrack->ptot() >= 0.2 && femtoTrack->ptot() < 3.4 &&
- TMath::Abs(femtoTrack->nSigmaPion()) < 3 &&
- femtoTrack->massSqr() >= -0.15 && femtoTrack->massSqr() < 0.1)
- {
- hdEdxQp[0][1]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
- hM2Qp[0][1]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
- hNsigQp[0][1]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
- hM2pt[0][1]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
- hNsigpt[0][1]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
- }
- // kaon
- if (femtoTrack->pt() >= 0.2 && femtoTrack->ptot() < 3.4 &&
- TMath::Abs(femtoTrack->nSigmaKaon()) < 3 &&
- femtoTrack->massSqr() >= 0.2 && femtoTrack->massSqr() < 0.32)
- {
- hdEdxQp[0][2]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
- hM2Qp[0][2]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
- hNsigQp[0][2]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
- hM2pt[0][2]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
- hNsigpt[0][2]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
- }
- // proton
- if (femtoTrack->ptot() >= 0.4 && femtoTrack->ptot() < 3.4 &&
- TMath::Abs(femtoTrack->nSigmaProton()) < 3 &&
- femtoTrack->massSqr() >= 0.75 && femtoTrack->massSqr() < 1.2)
- {
- hdEdxQp[0][3]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->dEdx());
- hM2Qp[0][3]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->massSqr());
- hNsigQp[0][3]->Fill(femtoTrack->charge() * femtoTrack->ptot(), femtoTrack->nSigmaPion());
- hM2pt[0][3]->Fill(femtoTrack->pt(), femtoTrack->massSqr());
- hNsigpt[0][3]->Fill(femtoTrack->pt(), femtoTrack->nSigmaPion());
- }
- } //if( isTofTrack() )
- } //for(Int_t iTrk=0; iTrk<nTracks; iTrk++)
- //Getting PsiEP with recentering + shift correction
- Psi_ReCenter=0.;
- Double_t mCos, mSin;
- if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
- {
- Psi2FullEP = 0.;
- Psi3FullEP = 0.;
- //FULL 2-nd order
- Psi_ReCenter = TMath::ATan2(fQv2FullEP.Y(), fQv2FullEP.X()) / 2.;
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax; iK++)
- {
- mCos = p_Cos2_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- mSin = p_Cos2_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
- }
- Psi2FullEP = GetEPAngle(Psi_ReCenter + dPsi, 2.);
- //FULL 3-nd order
- Psi_ReCenter = TMath::ATan2(fQv3FullEP.Y(), fQv3FullEP.X()) / 3.;
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax; iK++)
- {
- mCos = p_Cos3_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- mSin = p_Cos3_FULL_EP[VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
- }
- Psi3FullEP = GetEPAngle(Psi_ReCenter + dPsi, 3.);
- }
- for (int iGap=0; iGap<NEtaGaps; iGap++)
- {
- if (fQcountEast[iGap] > cutNcountQvGap &&
- fQcountWest[iGap] > cutNcountQvGap)
- {
- Psi2EastEP[iGap]=0.;
- Psi3EastEP[iGap]=0.;
- Psi2WestEP[iGap]=0.;
- Psi3WestEP[iGap]=0.;
- Psi2EastSP[iGap]=0.;
- Psi3EastSP[iGap]=0.;
- Psi2WestSP[iGap]=0.;
- Psi3WestSP[iGap]=0.;
- //EAST 2-nd order
- Psi_ReCenter = TMath::ATan2(fQv2EastEP[iGap].Y(),fQv2EastEP[iGap].X())/2.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- mCos = p_Cos2_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos2_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
- }
- Psi2EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
- Psi_ReCenter = TMath::ATan2(fQv2EastSP[iGap].Y(), fQv2EastSP[iGap].X()) / 2.;
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax; iK++)
- {
- mCos = p_Cos2_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos2_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
- }
- Psi2EastSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
- fQv2EastSP_Shift[iGap].Set(fQv2EastSP[iGap].Mod()*TMath::Cos(2.*Psi2EastSP[iGap]),fQv2EastSP[iGap].Mod()*TMath::Sin(2.*Psi2EastSP[iGap]));
- //WEST 2-nd order
- Psi_ReCenter = TMath::ATan2(fQv2WestEP[iGap].Y(),fQv2WestEP[iGap].X())/2.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- mCos = p_Cos2_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos2_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1./2.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff2.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff2.at(iK)*Psi_ReCenter));
- }
- Psi2WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,2.);
- Psi_ReCenter = TMath::ATan2(fQv2WestSP[iGap].Y(), fQv2WestSP[iGap].X()) / 2.;
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax; iK++)
- {
- mCos = p_Cos2_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos2_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
- }
- Psi2WestSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 2.);
- fQv2WestSP_Shift[iGap].Set(fQv2WestSP[iGap].Mod()*TMath::Cos(2.*Psi2WestSP[iGap]),fQv2WestSP[iGap].Mod()*TMath::Sin(2.*Psi2WestSP[iGap]));
- //EAST 3-nd order
- Psi_ReCenter = TMath::ATan2(fQv3EastEP[iGap].Y(),fQv3EastEP[iGap].X())/3.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- mCos = p_Cos3_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos3_EAST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
- }
- Psi3EastEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
- Psi_ReCenter = TMath::ATan2(fQv3EastSP[iGap].Y(), fQv3EastSP[iGap].X()) / 3.;
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax; iK++)
- {
- mCos = p_Cos3_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos3_EAST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
- }
- Psi3EastSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
- fQv3EastSP_Shift[iGap].Set(fQv3EastSP[iGap].Mod()*TMath::Cos(3.*Psi3EastSP[iGap]),fQv3EastSP[iGap].Mod()*TMath::Sin(3.*Psi3EastSP[iGap]));
- //WEST 3-nd order
- Psi_ReCenter = TMath::ATan2(fQv3WestEP[iGap].Y(),fQv3WestEP[iGap].X())/3.;
- dPsi = 0.;
- for (int iK=0; iK<NShiftOrderMax; iK++)
- {
- mCos = p_Cos3_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos3_WEST_EP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1./3.)*(2./(iK+1))*(-1.*mSin*
- TMath::Cos(shiftOrderCoeff3.at(iK)*Psi_ReCenter)+
- mCos*TMath::Sin(shiftOrderCoeff3.at(iK)*Psi_ReCenter));
- }
- Psi3WestEP[iGap] = AngleShift(Psi_ReCenter+dPsi,3.);
- Psi_ReCenter = TMath::ATan2(fQv3WestSP[iGap].Y(), fQv3WestSP[iGap].X()) / 3.;
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax; iK++)
- {
- mCos = p_Cos3_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- mSin = p_Cos3_WEST_SP[VtxSign][iK][(Int_t)event->cent16()][iGap][(Double_t)event->runId()];
- dPsi += (1. / 3.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff3.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff3.at(iK) * Psi_ReCenter));
- }
- Psi3WestSP[iGap] = AngleShift(Psi_ReCenter + dPsi, 3.);
- fQv3WestSP_Shift[iGap].Set(fQv3WestSP[iGap].Mod()*TMath::Cos(3.*Psi3WestSP[iGap]),fQv3WestSP[iGap].Mod()*TMath::Sin(3.*Psi3WestSP[iGap]));
- }
- }
- // 2nd harmonics TPC
- if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
- {
- hTPCFullQxRaw[0]->Fill(fQv2FullEPRaw.X(),event->cent16());
- hTPCFullQyRaw[0]->Fill(fQv2FullEPRaw.X(),event->cent16());
- hTPCFullPsiRaw[0]->Fill(GetEPAngle(fQv2FullEPRaw,2), event->cent16());
- }
- if (fQcountEast[0] > cutNcountQvGap &&
- fQcountWest[0] > cutNcountQvGap)
- {
- hTPCEastQxRaw[0]->Fill(fQv2EastEPRaw[0].X(),event->cent16());
- hTPCEastQyRaw[0]->Fill(fQv2EastEPRaw[0].Y(),event->cent16());
- hTPCWestQxRaw[0]->Fill(fQv2WestEPRaw[0].X(),event->cent16());
- hTPCWestQyRaw[0]->Fill(fQv2WestEPRaw[0].Y(),event->cent16());
- hTPCEastPsiRaw[0]->Fill(GetEPAngle(fQv2EastEPRaw[0],2), event->cent16());
- hTPCWestPsiRaw[0]->Fill(GetEPAngle(fQv2WestEPRaw[0],2), event->cent16());
- }
- if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
- {
- hTPCFullQxRec[0]->Fill(fQv2FullEP.X(),event->cent16());
- hTPCFullQyRec[0]->Fill(fQv2FullEP.X(),event->cent16());
- hTPCFullPsiRec[0]->Fill(GetEPAngle(fQv2FullEP,2), event->cent16());
- }
- if (fQcountEast[0] > cutNcountQvGap &&
- fQcountWest[0] > cutNcountQvGap)
- {
- hTPCEastQxRec[0]->Fill(fQv2EastEP[0].X(),event->cent16());
- hTPCEastQyRec[0]->Fill(fQv2EastEP[0].Y(),event->cent16());
- hTPCWestQxRec[0]->Fill(fQv2WestEP[0].X(),event->cent16());
- hTPCWestQyRec[0]->Fill(fQv2WestEP[0].Y(),event->cent16());
- hTPCEastPsiRec[0]->Fill(GetEPAngle(fQv2EastEP[0],2), event->cent16());
- hTPCWestPsiRec[0]->Fill(GetEPAngle(fQv2WestEP[0],2), event->cent16());
- }
- if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
- {
- hTPCFullPsiShift[0]->Fill(GetEPAngle(Psi2FullEP,2), event->cent16());
- // East
- pQx2TPCvsRun[0]->Fill(event->runId(),fQv2EastEP[0].X());
- pQy2TPCvsRun[0]->Fill(event->runId(),fQv2EastEP[0].Y());
- pQx3TPCvsRun[0]->Fill(event->runId(),fQv3EastEP[0].X());
- pQy3TPCvsRun[0]->Fill(event->runId(),fQv3EastEP[0].Y());
- pCos2TPCvsRun[0][0]->Fill(event->runId(),TMath::Cos(1*(Psi2EastEP[0])));
- pCos2TPCvsRun[1][0]->Fill(event->runId(),TMath::Cos(2*(Psi2EastEP[0])));
- pCos2TPCvsRun[2][0]->Fill(event->runId(),TMath::Cos(3*(Psi2EastEP[0])));
- pCos2TPCvsRun[3][0]->Fill(event->runId(),TMath::Cos(4*(Psi2EastEP[0])));
- pCos3TPCvsRun[0][0]->Fill(event->runId(),TMath::Cos(1*(Psi3EastEP[0])));
- pCos3TPCvsRun[1][0]->Fill(event->runId(),TMath::Cos(2*(Psi3EastEP[0])));
- pCos3TPCvsRun[2][0]->Fill(event->runId(),TMath::Cos(3*(Psi3EastEP[0])));
- pCos3TPCvsRun[3][0]->Fill(event->runId(),TMath::Cos(4*(Psi3EastEP[0])));
- pSin2TPCvsRun[0][0]->Fill(event->runId(),TMath::Sin(1*(Psi2EastEP[0])));
- pSin2TPCvsRun[1][0]->Fill(event->runId(),TMath::Sin(2*(Psi2EastEP[0])));
- pSin2TPCvsRun[2][0]->Fill(event->runId(),TMath::Sin(3*(Psi2EastEP[0])));
- pSin2TPCvsRun[3][0]->Fill(event->runId(),TMath::Sin(4*(Psi2EastEP[0])));
- pSin3TPCvsRun[0][0]->Fill(event->runId(),TMath::Sin(1*(Psi3EastEP[0])));
- pSin3TPCvsRun[1][0]->Fill(event->runId(),TMath::Sin(2*(Psi3EastEP[0])));
- pSin3TPCvsRun[2][0]->Fill(event->runId(),TMath::Sin(3*(Psi3EastEP[0])));
- pSin3TPCvsRun[3][0]->Fill(event->runId(),TMath::Sin(4*(Psi3EastEP[0])));
- // West
- pQx2TPCvsRun[1]->Fill(event->runId(),fQv2WestEP[0].X());
- pQy2TPCvsRun[1]->Fill(event->runId(),fQv2WestEP[0].Y());
- pQx3TPCvsRun[1]->Fill(event->runId(),fQv3WestEP[0].X());
- pQy3TPCvsRun[1]->Fill(event->runId(),fQv3WestEP[0].Y());
- pCos2TPCvsRun[0][1]->Fill(event->runId(),TMath::Cos(1*(Psi2WestEP[0])));
- pCos2TPCvsRun[1][1]->Fill(event->runId(),TMath::Cos(2*(Psi2WestEP[0])));
- pCos2TPCvsRun[2][1]->Fill(event->runId(),TMath::Cos(3*(Psi2WestEP[0])));
- pCos2TPCvsRun[3][1]->Fill(event->runId(),TMath::Cos(4*(Psi2WestEP[0])));
- pCos3TPCvsRun[0][1]->Fill(event->runId(),TMath::Cos(1*(Psi3WestEP[0])));
- pCos3TPCvsRun[1][1]->Fill(event->runId(),TMath::Cos(2*(Psi3WestEP[0])));
- pCos3TPCvsRun[2][1]->Fill(event->runId(),TMath::Cos(3*(Psi3WestEP[0])));
- pCos3TPCvsRun[3][1]->Fill(event->runId(),TMath::Cos(4*(Psi3WestEP[0])));
- pSin2TPCvsRun[0][1]->Fill(event->runId(),TMath::Sin(1*(Psi2WestEP[0])));
- pSin2TPCvsRun[1][1]->Fill(event->runId(),TMath::Sin(2*(Psi2WestEP[0])));
- pSin2TPCvsRun[2][1]->Fill(event->runId(),TMath::Sin(3*(Psi2WestEP[0])));
- pSin2TPCvsRun[3][1]->Fill(event->runId(),TMath::Sin(4*(Psi2WestEP[0])));
- pSin3TPCvsRun[0][1]->Fill(event->runId(),TMath::Sin(1*(Psi3WestEP[0])));
- pSin3TPCvsRun[1][1]->Fill(event->runId(),TMath::Sin(2*(Psi3WestEP[0])));
- pSin3TPCvsRun[2][1]->Fill(event->runId(),TMath::Sin(3*(Psi3WestEP[0])));
- pSin3TPCvsRun[3][1]->Fill(event->runId(),TMath::Sin(4*(Psi3WestEP[0])));
- }
- if (fQcountEast[0] > cutNcountQvGap &&
- fQcountWest[0] > cutNcountQvGap)
- {
- hTPCEastPsiShift[0]->Fill(GetEPAngle(Psi2EastEP[0],2), event->cent16());
- hTPCWestPsiShift[0]->Fill(GetEPAngle(Psi2WestEP[0],2), event->cent16());
- }
- // 3rd harmonics TPC
- if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
- {
- hTPCFullQxRaw[1]->Fill(fQv3FullEPRaw.X(),event->cent16());
- hTPCFullQyRaw[1]->Fill(fQv3FullEPRaw.X(),event->cent16());
- hTPCFullPsiRaw[1]->Fill(GetEPAngle(fQv3FullEPRaw,3), event->cent16());
- }
- if (fQcountEast[0] > cutNcountQvGap &&
- fQcountWest[0] > cutNcountQvGap)
- {
- hTPCEastQxRaw[1]->Fill(fQv3EastEPRaw[0].X(),event->cent16());
- hTPCEastQyRaw[1]->Fill(fQv3EastEPRaw[0].Y(),event->cent16());
- hTPCWestQxRaw[1]->Fill(fQv3WestEPRaw[0].X(),event->cent16());
- hTPCWestQyRaw[1]->Fill(fQv3WestEPRaw[0].Y(),event->cent16());
- hTPCEastPsiRaw[1]->Fill(GetEPAngle(fQv3EastEPRaw[0],3), event->cent16());
- hTPCWestPsiRaw[1]->Fill(GetEPAngle(fQv3WestEPRaw[0],3), event->cent16());
- }
- if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
- {
- hTPCFullQxRec[1]->Fill(fQv3FullEP.X(),event->cent16());
- hTPCFullQyRec[1]->Fill(fQv3FullEP.X(),event->cent16());
- hTPCFullPsiRec[1]->Fill(GetEPAngle(fQv3FullEP,3), event->cent16());
- }
- if (fQcountEast[0] > cutNcountQvGap &&
- fQcountWest[0] > cutNcountQvGap)
- {
- hTPCEastQxRec[1]->Fill(fQv3EastEP[0].X(),event->cent16());
- hTPCEastQyRec[1]->Fill(fQv3EastEP[0].Y(),event->cent16());
- hTPCWestQxRec[1]->Fill(fQv3WestEP[0].X(),event->cent16());
- hTPCWestQyRec[1]->Fill(fQv3WestEP[0].Y(),event->cent16());
- hTPCEastPsiRec[1]->Fill(GetEPAngle(fQv3EastEP[0],3), event->cent16());
- hTPCWestPsiRec[1]->Fill(GetEPAngle(fQv3WestEP[0],3), event->cent16());
- }
- if (fQcountFullWest + fQcountFullEast > cutNcountQvFull)
- {
- hTPCFullPsiShift[1]->Fill(GetEPAngle(Psi3FullEP,3), event->cent16());
- }
- if (fQcountEast[0] > cutNcountQvGap &&
- fQcountWest[0] > cutNcountQvGap)
- {
- hTPCEastPsiShift[1]->Fill(GetEPAngle(Psi3EastEP[0],3), event->cent16());
- hTPCWestPsiShift[1]->Fill(GetEPAngle(Psi3WestEP[0],3), event->cent16());
- }
- } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
- femtoReader->Finish();
- TFile *output = new TFile(outFile, "recreate");
- for (int iHarm=0; iHarm<fNharmonics+1; iHarm++)
- {
- // Raw
- hBBCFullQxRaw[iHarm]->Write();
- hBBCFullQyRaw[iHarm]->Write();
- hBBCEastQxRaw[iHarm]->Write();
- hBBCEastQyRaw[iHarm]->Write();
- hBBCWestQxRaw[iHarm]->Write();
- hBBCWestQyRaw[iHarm]->Write();
- hBBCFullPsiRaw[iHarm]->Write();
- hBBCEastPsiRaw[iHarm]->Write();
- hBBCWestPsiRaw[iHarm]->Write();
- // w/ Gain correction
- hBBCFullQxGain[iHarm]->Write();
- hBBCFullQyGain[iHarm]->Write();
- hBBCEastQxGain[iHarm]->Write();
- hBBCEastQyGain[iHarm]->Write();
- hBBCWestQxGain[iHarm]->Write();
- hBBCWestQyGain[iHarm]->Write();
- hBBCFullPsiGain[iHarm]->Write();
- hBBCEastPsiGain[iHarm]->Write();
- hBBCWestPsiGain[iHarm]->Write();
- // w/ Gain + recentering
- hBBCFullQxRec[iHarm]->Write();
- hBBCFullQyRec[iHarm]->Write();
- hBBCEastQxRec[iHarm]->Write();
- hBBCEastQyRec[iHarm]->Write();
- hBBCWestQxRec[iHarm]->Write();
- hBBCWestQyRec[iHarm]->Write();
- hBBCFullPsiRec[iHarm]->Write();
- hBBCEastPsiRec[iHarm]->Write();
- hBBCWestPsiRec[iHarm]->Write();
- // w/ Gain + recentering + shift
- hBBCFullPsiShift[iHarm]->Write();
- hBBCEastPsiShift[iHarm]->Write();
- hBBCWestPsiShift[iHarm]->Write();
- }
- for (int iHarm=0; iHarm<fNharmonics; iHarm++)
- {
- // Raw
- hTPCFullQxRaw[iHarm]->Write();
- hTPCFullQyRaw[iHarm]->Write();
- hTPCEastQxRaw[iHarm]->Write();
- hTPCEastQyRaw[iHarm]->Write();
- hTPCWestQxRaw[iHarm]->Write();
- hTPCWestQyRaw[iHarm]->Write();
- hTPCFullPsiRaw[iHarm]->Write();
- hTPCEastPsiRaw[iHarm]->Write();
- hTPCWestPsiRaw[iHarm]->Write();
- // w/ recentering
- hTPCFullQxRec[iHarm]->Write();
- hTPCFullQyRec[iHarm]->Write();
- hTPCEastQxRec[iHarm]->Write();
- hTPCEastQyRec[iHarm]->Write();
- hTPCWestQxRec[iHarm]->Write();
- hTPCWestQyRec[iHarm]->Write();
- hTPCFullPsiRec[iHarm]->Write();
- hTPCEastPsiRec[iHarm]->Write();
- hTPCWestPsiRec[iHarm]->Write();
- // w/ recentering + shift
- hTPCFullPsiShift[iHarm]->Write();
- hTPCEastPsiShift[iHarm]->Write();
- hTPCWestPsiShift[iHarm]->Write();
- }
- for (int i = 0; i < 2; i++)
- {
- //Eventwise
- hVtxXY[i]->Write();
- hVtxZ[i]->Write();
- hVtxVpdZ[i]->Write();
- hVpdZ[i]->Write();
- pVtxXvsRun[i]->Write();
- pVtxYvsRun[i]->Write();
- pVtxZvsRun[i]->Write();
- pVpdZvsRun[i]->Write();
- pVtxZVpdZvsRun[i]->Write();
- pRefMultvsRun[i]->Write();
- pTofMultvsRun[i]->Write();
- pTofMatchvsRun[i]->Write();
- pAdcvsRun[i]->Write();
- pNtrPrimvsRun[i]->Write();
- pNtrGlobvsRun[i]->Write();
-
- //Trackwise
- hNhitsFit[i]->Write();
- hNhitsPoss[i]->Write();
- hNhitsRatio[i]->Write();
- hPt[i]->Write();
- hPhi[i]->Write();
- hPtot[i]->Write();
- hEta[i]->Write();
- hDCA[i]->Write();
- pPtvsRun[i]->Write();
- pEtavsRun[i]->Write();
- pPhivsRun[i]->Write();
- pNhitsFitvsRun[i]->Write();
- pNhitsPossvsRun[i]->Write();
- pDCAvsRun[i]->Write();
- pdEdxvsRun[i]->Write();
- pQx2TPCvsRun[i]->Write();
- pQy2TPCvsRun[i]->Write();
- pQx3TPCvsRun[i]->Write();
- pQy3TPCvsRun[i]->Write();
- pQx2BBCvsRun[i]->Write();
- pQy2BBCvsRun[i]->Write();
- for (int iCent=0;iCent<2*fNharmonics;iCent++)
- {
- pCos2TPCvsRun[iCent][i]->Write();
- pSin2TPCvsRun[iCent][i]->Write();
- pCos3TPCvsRun[iCent][i]->Write();
- pSin3TPCvsRun[iCent][i]->Write();
- pCosBBCvsRun[iCent][i]->Write();
- pSinBBCvsRun[iCent][i]->Write();
- }
- //PID
- hnSigM2_lowpt[i]->Write();
- hnSigM2_highpt[i]->Write();
- for (int iPID = 0; iPID < 4; iPID++)
- {
- hdEdxQp[i][iPID]->Write();
- hM2Qp[i][iPID]->Write();
- hNsigQp[i][iPID]->Write();
- hM2pt[i][iPID]->Write();
- hNsigpt[i][iPID]->Write();
- }
- }
- for (int iPID = 0; iPID < 4; iPID++)
- {
- for (int iCent=0;iCent<16;iCent++)
- {
- h_pt_spectra[iPID][iCent]->Write();
- }
- }
- hNevVsRun->Write();
- output->Close();
- std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
- << std::endl;
- timer.Stop();
- timer.Print();
- }
- Bool_t isGoodEvent(StFemtoEvent *const &event)
- {
- if (!event)
- return false;
- if (event == nullptr)
- return false;
- if (event->primaryVertex().Perp() > cutVtxR)
- return false;
- if (TMath::Abs(event->primaryVertex().Z()) > cutVtxZEnergy.at(energy))
- return false;
- if ((energy == 200.) && TMath::Abs(event->primaryVertex().Z() - event->vpdVz()) > cutVpdVz)
- return false;
- return true;
- }
- Bool_t isGoodTrack(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
- {
- if (!track)
- return false;
- // if (!track->isPrimary()) return false;
- if (track->nHitsFit() < cutNhits)
- return false;
- if (track->nHitsPoss() <= cutNhitsPoss)
- return false;
- if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
- return false;
- if (TMath::Abs(track->eta()) >= cutEta)
- return false;
- if (track->pt() <= cutPtMin.at(_energy))
- return false;
- if (track->pt() > cutPtMax)
- return false;
- if (track->ptot() > cutPMax)
- return false;
- if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
- return false;
- return true;
- }
- Bool_t isGoodTrackFlow(StFemtoTrack *const &track, Float_t _energy, TVector3 pVtx)
- {
- if (!track)
- return false;
- // if (!track->isPrimary()) return false;
- if (track->nHitsFit() < cutNhits)
- return false;
- if (track->nHitsPoss() <= cutNhitsPoss)
- return false;
- if ((Double_t)track->nHitsFit() / track->nHitsPoss() < cutNhitsRatio)
- return false;
- if (TMath::Abs(track->eta()) >= cutEta)
- return false;
- if (track->pt() <= cutPtMin.at(_energy))
- return false;
- //if (track->pt() > cutPtMax) return false;
- if (track->ptot() > cutPMax)
- return false;
- if (track->gDCA(pVtx).Mag() >= cutDCA.at(_energy))
- return false;
- return true;
- }
- Double_t GetWeight(StFemtoTrack *const &track)
- {
- Double_t weight;
- if (track->pt() <= cutPtWeightEP)
- {
- weight = track->pt();
- }
- else
- {
- weight = cutPtWeightEP;
- }
- return weight;
- }
- TVector2 CalcQvector(StFemtoTrack *const &track, Int_t _harm)
- {
- TVector2 qv(0., 0.);
- qv.Set(TMath::Cos(_harm * track->phi()), TMath::Sin(_harm * track->phi()));
- return qv;
- }
- Double_t AngleShift(Double_t Psi, Double_t order)
- {
- Double_t PsiCorr = Psi;
- if (PsiCorr > TMath::Pi() / order)
- {
- PsiCorr = Psi - 2. * TMath::Pi() / order;
- }
- if (PsiCorr < -1. * TMath::Pi() / order)
- {
- PsiCorr = Psi + 2. * TMath::Pi() / order;
- }
- return PsiCorr;
- }
- Int_t GetVzBin(Double_t vtxZ)
- {
- for (Int_t i = 0; i < fArraySize; i++)
- {
- if (vtxZ >= VtxBins[i] && vtxZ < VtxBins[i + 1])
- return i;
- }
- return -1;
- }
- //--------------------------------------------------------------------------------------------------------------------------//
- //this function simply connects the gain values read in to the BBC azimuthal distribution
- //since tiles 7 and 9 (+ 13 and 15) share a gain value it is ambiguous how to assign the geometry here
- //I prefer assigning the angle between the tiles thus "greying out" the adcs.
- //Others have assigned all of the adc to one (exclusive) or the the other.
- Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId)
- {
- //float GetPhiInBBC(int eastWest, int bbcN) { //tileId=0 to 23
- const float Pi = TMath::Pi();
- const float Phi_div = Pi / 6;
- float bbc_phi = Phi_div;
- switch(tileId) {
- case 0: bbc_phi = 3.*Phi_div;
- break;
- case 1: bbc_phi = Phi_div;
- break;
- case 2: bbc_phi = -1.*Phi_div;
- break;
- case 3: bbc_phi = -3.*Phi_div;
- break;
- case 4: bbc_phi = -5.*Phi_div;
- break;
- case 5: bbc_phi = 5.*Phi_div;
- break;
- //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
- case 6: bbc_phi = 3.*Phi_div;
- break;
- case 7: bbc_phi = 3.*Phi_div;
- break;
- case 8: bbc_phi = Phi_div;
- break;
- case 9: bbc_phi = 0.;
- break;
- case 10: bbc_phi = -1.*Phi_div;
- break;
- //case 11: bbc_phi = (mRndm.Rndm() > 0.5) ? -2.*Phi_div:-4.*Phi_div; //tiles 13 and 15 are gained together
- case 11: bbc_phi = -3.*Phi_div;
- break;
- case 12: bbc_phi = -3.*Phi_div;
- break;
- case 13: bbc_phi = -5.*Phi_div;
- break;
- case 14: bbc_phi = Pi;
- break;
- case 15: bbc_phi = 5.*Phi_div;
- break;
- }
- //if we're looking at the east BBC we need to flip around x in the STAR coordinates,
- //a line parallel to the beam would go through tile 1 on the W BBC and tile 3 on the
- if(0 == eastWest){
- if (bbc_phi > -0.001){ //this is not a >= since we are talking about finite adcs -- not to important
- bbc_phi = Pi - bbc_phi;
- }
- else {
- bbc_phi= -Pi - bbc_phi;
- }
- }
- if(bbc_phi < 0.0) bbc_phi += 2.*Pi;
- if(bbc_phi > 2.*Pi) bbc_phi -= 2.*Pi;
- return bbc_phi;
- }
- Double_t GetEPAngle(TVector2 qvect, Double_t order)
- {
- Double_t angle = -9999.;
- if(qvect.Mod())
- {
- angle = qvect.Phi();
- if(angle < 0.0) {angle += 2.*TMath::Pi()/(order);}
- if(angle > 2.*TMath::Pi()/(order)) {angle -= 2.*TMath::Pi()/(order);}
- }
- return angle;
- }
- Double_t GetEPAngle(Double_t angle, Double_t order)
- {
- if(angle < 0.0) {angle += 2.*TMath::Pi()/(order);}
- if(angle > 2.*TMath::Pi()/(order)) {angle -= 2.*TMath::Pi()/(order);}
- return angle;
- }
- Bool_t isGoodPID(StFemtoTrack *const &track)
- {
- if (!track->isTofTrack()) return false;
- if (track->massSqr() < cutMass2Min) return false;
- //NhitsDedx cut applied in StFemtoDst?
- // ToFYLocal cut applied in StFemtoDst?
- return true;
- }
|