1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071727374757677787980818283848586878889909192939495969798991001011021031041051061071081091101111121131141151161171181191201211221231241251261271281291301311321331341351361371381391401411421431441451461471481491501511521531541551561571581591601611621631641651661671681691701711721731741751761771781791801811821831841851861871881891901911921931941951961971981992002012022032042052062072082092102112122132142152162172182192202212222232242252262272282292302312322332342352362372382392402412422432442452462472482492502512522532542552562572582592602612622632642652662672682692702712722732742752762772782792802812822832842852862872882892902912922932942952962972982993003013023033043053063073083093103113123133143153163173183193203213223233243253263273283293303313323333343353363373383393403413423433443453463473483493503513523533543553563573583593603613623633643653663673683693703713723733743753763773783793803813823833843853863873883893903913923933943953963973983994004014024034044054064074084094104114124134144154164174184194204214224234244254264274284294304314324334344354364374384394404414424434444454464474484494504514524534544554564574584594604614624634644654664674684694704714724734744754764774784794804814824834844854864874884894904914924934944954964974984995005015025035045055065075085095105115125135145155165175185195205215225235245255265275285295305315325335345355365375385395405415425435445455465475485495505515525535545555565575585595605615625635645655665675685695705715725735745755765775785795805815825835845855865875885895905915925935945955965975985996006016026036046056066076086096106116126136146156166176186196206216226236246256266276286296306316326336346356366376386396406416426436446456466476486496506516526536546556566576586596606616626636646656666676686696706716726736746756766776786796806816826836846856866876886896906916926936946956966976986997007017027037047057067077087097107117127137147157167177187197207217227237247257267277287297307317327337347357367377387397407417427437447457467477487497507517527537547557567577587597607617627637647657667677687697707717727737747757767777787797807817827837847857867877887897907917927937947957967977987998008018028038048058068078088098108118128138148158168178188198208218228238248258268278288298308318328338348358368378388398408418428438448458468478488498508518528538548558568578588598608618628638648658668678688698708718728738748758768778788798808818828838848858868878888898908918928938948958968978988999009019029039049059069079089099109119129139149159169179189199209219229239249259269279289299309319329339349359369379389399409419429439449459469479489499509519529539549559569579589599609619629639649659669679689699709719729739749759769779789799809819829839849859869879889899909919929939949959969979989991000100110021003100410051006100710081009101010111012101310141015101610171018101910201021102210231024102510261027102810291030103110321033103410351036103710381039104010411042104310441045104610471048104910501051105210531054105510561057105810591060106110621063106410651066106710681069107010711072107310741075107610771078107910801081108210831084108510861087108810891090109110921093109410951096109710981099110011011102110311041105110611071108110911101111111211131114111511161117111811191120112111221123112411251126112711281129113011311132113311341135113611371138113911401141114211431144114511461147114811491150115111521153115411551156115711581159116011611162116311641165116611671168116911701171117211731174117511761177117811791180118111821183118411851186118711881189119011911192119311941195119611971198119912001201120212031204120512061207120812091210121112121213121412151216121712181219122012211222122312241225122612271228122912301231123212331234123512361237123812391240124112421243124412451246124712481249125012511252125312541255125612571258125912601261126212631264126512661267126812691270127112721273127412751276127712781279128012811282128312841285128612871288128912901291129212931294129512961297129812991300130113021303130413051306130713081309131013111312131313141315131613171318131913201321132213231324132513261327132813291330133113321333133413351336133713381339134013411342134313441345134613471348134913501351135213531354135513561357135813591360136113621363136413651366136713681369137013711372137313741375137613771378137913801381138213831384138513861387138813891390139113921393139413951396139713981399140014011402140314041405140614071408140914101411141214131414141514161417141814191420142114221423142414251426142714281429143014311432143314341435143614371438143914401441144214431444144514461447144814491450145114521453145414551456145714581459146014611462146314641465146614671468146914701471147214731474147514761477147814791480148114821483148414851486148714881489149014911492149314941495149614971498149915001501150215031504150515061507150815091510151115121513151415151516151715181519152015211522152315241525152615271528152915301531153215331534153515361537153815391540154115421543154415451546154715481549155015511552155315541555155615571558155915601561156215631564156515661567156815691570157115721573157415751576157715781579158015811582158315841585158615871588158915901591159215931594159515961597159815991600160116021603160416051606160716081609161016111612161316141615161616171618161916201621162216231624162516261627162816291630163116321633163416351636163716381639164016411642164316441645164616471648164916501651165216531654165516561657165816591660166116621663166416651666166716681669167016711672167316741675167616771678167916801681168216831684168516861687168816891690169116921693169416951696169716981699170017011702170317041705170617071708170917101711171217131714171517161717171817191720172117221723172417251726172717281729173017311732173317341735173617371738173917401741174217431744174517461747174817491750175117521753175417551756175717581759176017611762176317641765176617671768176917701771177217731774177517761777177817791780178117821783178417851786178717881789179017911792179317941795179617971798179918001801180218031804180518061807180818091810181118121813181418151816181718181819182018211822182318241825182618271828182918301831183218331834183518361837183818391840184118421843184418451846184718481849185018511852185318541855185618571858185918601861186218631864186518661867186818691870187118721873187418751876187718781879188018811882188318841885188618871888188918901891189218931894189518961897189818991900190119021903190419051906190719081909191019111912191319141915191619171918191919201921192219231924192519261927192819291930193119321933193419351936193719381939194019411942194319441945194619471948194919501951195219531954195519561957195819591960196119621963196419651966196719681969197019711972197319741975197619771978197919801981198219831984198519861987198819891990199119921993199419951996199719981999200020012002200320042005200620072008200920102011201220132014201520162017201820192020202120222023202420252026202720282029203020312032203320342035203620372038203920402041204220432044204520462047204820492050205120522053205420552056205720582059206020612062206320642065206620672068206920702071207220732074207520762077207820792080208120822083208420852086208720882089209020912092209320942095209620972098209921002101210221032104210521062107210821092110211121122113211421152116211721182119212021212122212321242125212621272128212921302131213221332134213521362137213821392140214121422143214421452146214721482149215021512152215321542155215621572158215921602161216221632164216521662167216821692170217121722173217421752176217721782179218021812182218321842185218621872188218921902191219221932194219521962197219821992200220122022203220422052206220722082209221022112212221322142215221622172218221922202221222222232224222522262227222822292230223122322233223422352236223722382239224022412242224322442245224622472248224922502251225222532254225522562257225822592260226122622263226422652266226722682269227022712272227322742275227622772278227922802281228222832284228522862287228822892290229122922293229422952296229722982299230023012302230323042305230623072308230923102311231223132314231523162317231823192320232123222323232423252326232723282329233023312332233323342335233623372338233923402341234223432344234523462347234823492350235123522353235423552356235723582359236023612362236323642365236623672368236923702371237223732374237523762377237823792380238123822383238423852386238723882389239023912392239323942395239623972398239924002401240224032404240524062407240824092410241124122413241424152416241724182419242024212422242324242425242624272428242924302431243224332434243524362437243824392440244124422443244424452446244724482449245024512452245324542455245624572458245924602461246224632464246524662467246824692470247124722473247424752476247724782479248024812482248324842485248624872488248924902491249224932494249524962497249824992500250125022503250425052506250725082509251025112512251325142515251625172518251925202521252225232524252525262527252825292530 |
- /**
- * \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"
- // 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);
- Int_t GetVzBin(Double_t vtxZ);
- Float_t BBC_GetPhi(const Int_t eastWest, const Int_t tileId);
- Bool_t isBadRun(StFemtoEvent *const &event);
- Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip);
- //_________________
- void FemtoDstAnalyzer_FlowChargedHadrons(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 *resFile = "",
- const Char_t *gainBBC = "",
- const Char_t *reCentFileBBC = "",
- const Char_t *shiftFileBBC = "",
- const Char_t *reCentFileZDC = "",
- const Char_t *shiftFileZDC = "",
- const Char_t *resFitFile = "")
- {
- 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;
- }
- }
- }
- //Read resolution
- TFile *fiRes = new TFile(resFile, "read");
- TProfile *p_res2_EP[NEtaGaps];
- TProfile *p_res3_EP[NEtaGaps];
- TProfile *p_res2_SP[NEtaGaps];
- TProfile *p_res3_SP[NEtaGaps];
- // TProfile *p_resBBC_EP[fNharmonics+1];
- TProfile *p_res2_EP_3sub_AB;
- TProfile *p_res2_EP_3sub_AC;
- TProfile *p_res2_EP_3sub_BC;
- p_res2_EP_3sub_AB = (TProfile *)fiRes->Get("p_res2_EP_3sub_AB");
- p_res2_EP_3sub_AC = (TProfile *)fiRes->Get("p_res2_EP_3sub_AC");
- p_res2_EP_3sub_BC = (TProfile *)fiRes->Get("p_res2_EP_3sub_BC");
- for (int iGap = 0; iGap < NEtaGaps; iGap++)
- {
- p_res2_EP[iGap] = (TProfile *)fiRes->Get(Form("p_res2_EP_gap%i", iGap));
- p_res3_EP[iGap] = (TProfile *)fiRes->Get(Form("p_res3_EP_gap%i", iGap));
- p_res2_SP[iGap] = (TProfile *)fiRes->Get(Form("p_res2_SP_gap%i", iGap));
- p_res3_SP[iGap] = (TProfile *)fiRes->Get(Form("p_res3_SP_gap%i", iGap));
- }
- // for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- // {
- // p_resBBC_EP[iHarm] = (TProfile*) fiRes->Get(Form("p_resBBC_EP%i", iHarm));
- // }
- //Read resolution fit
- TFile *fiResFit = new TFile(resFitFile, "read");
- TH1D *hResBBCFULL[fNharmonics + 1];
- TH1D *hResBBCHalf[fNharmonics + 1];
- TH1D *hResZDCFULL[fNharmonics + 1];
- TH1D *hResZDCHalf[fNharmonics + 1];
- for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- {
- hResBBCFULL[iHarm] = (TH1D *)fiResFit->Get(Form("hResBBCFULL%i", iHarm));
- hResBBCHalf[iHarm] = (TH1D *)fiResFit->Get(Form("hResBBCHalf%i", iHarm));
- }
- for (int iHarm = 0; iHarm < fNharmonics; iHarm++)
- {
- hResZDCFULL[iHarm] = (TH1D *)fiResFit->Get(Form("hResZDCFULL%i", iHarm));
- hResZDCHalf[iHarm] = (TH1D *)fiResFit->Get(Form("hResZDCHalf%i", iHarm));
- }
- //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;
- }
- }
- }
- //Manage recentering ZDC
- std::map<Double_t, Double_t> p_ZDC_Qx_Full_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Qy_Full_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Qx_East_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Qy_East_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Qx_West_EP[fNharmonics][fArraySize][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Qy_West_EP[fNharmonics][fArraySize][fNCent];
- TFile *fiReCentZDC = new TFile(reCentFileZDC, "read");
- std::cout << "Reading file for ZDC 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 *)fiReCentZDC->Get(Form("p_ZDC_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_ZDC_Qx_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_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_ZDC_Qy_Full_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_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_ZDC_Qx_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_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_ZDC_Qy_East_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_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_ZDC_Qx_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiReCentZDC->Get(Form("p_ZDC_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_ZDC_Qy_West_EP[iHarm][iVtx][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- }
- }
- //Manage ShiftCorr ZDC
- std::map<Double_t, Double_t> p_ZDC_Cos_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Sin_Full_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Cos_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Sin_East_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Cos_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- std::map<Double_t, Double_t> p_ZDC_Sin_West_EP[fNharmonics + 1][fArraySize][NShiftOrderMax * 3][fNCent];
- TFile *fiShiftZDC = new TFile(shiftFileZDC, "read");
- fiShiftZDC->cd();
- std::cout << "Reading file for ZDC 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 *)fiShiftZDC->Get(Form("p_ZDC_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_ZDC_Cos_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_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_ZDC_Sin_Full_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_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_ZDC_Cos_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_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_ZDC_Sin_East_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_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_ZDC_Cos_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- prof = (TProfile2D *)fiShiftZDC->Get(Form("p_ZDC_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_ZDC_Sin_West_EP[iHarm][iVtx][iShift][iCent][(Double_t)(runId_min + i)] = binCont;
- }
- }
- delete prof;
- }
- }
- }
- //Manage output TProfiles
- //TProfile *p_v2_EP[9][NEtaGaps];
- //TProfile *p_v3_EP[9][NEtaGaps];
- //TProfile *p_v2cent_EP[NEtaGaps];
- //TProfile *p_v3cent_EP[NEtaGaps];
- TProfile2D *p_v2cent_EP_ptgaps[NEtaGaps];
- TProfile2D *p_v3cent_EP_ptgaps[NEtaGaps];
- TProfile2D *p_v2cent_SP_ptgaps[NEtaGaps];
- TProfile2D *p_v3cent_SP_ptgaps[NEtaGaps];
- TProfile2D *p_v2WEST_EP_ptgaps[NEtaGaps];
- TProfile2D *p_v3WEST_EP_ptgaps[NEtaGaps];
- TProfile2D *p_v2EAST_EP_ptgaps[NEtaGaps];
- TProfile2D *p_v3EAST_EP_ptgaps[NEtaGaps];
- TProfile2D *p_vnBBC[fNharmonics + 1];
- // TProfile2D *p_vnZDC[fNharmonics+1];
- // TProfile2D *p_vnBBCEast[fNharmonics+1];
- // TProfile2D *p_vnBBCWest[fNharmonics+1];
- TProfile2D *p_v2run_EP[NEtaGaps][16]; // runId vs pt
- TProfile2D *p_v3run_EP[NEtaGaps][16]; // runId vs pt
- TProfile2D *p_v2run_SP[NEtaGaps][16]; // runId vs pt
- TProfile2D *p_v3run_SP[NEtaGaps][16]; // runId vs pt
- TProfile2D *p_v2run_BBC[16]; // runId vs pt
- TProfile2D *p_v2run_ZDC[16]; // runId vs pt
- TProfile2D *p_v2run_ZDC_East[16]; // runId vs pt
- TProfile2D *p_v2run_ZDC_West[16]; // runId vs pt
- //TH2D *p_ptrun[16]; // to calculate <pT> for each pT bin
- //TH2D *p1_ptrun[16]; // to calculate <pT> for each pT bin
- /*TProfile2D *p_v2cent_EP_3subTPC;
- TProfile2D *p_v2cent_EP_3subBBCEast;
- TProfile2D *p_v2cent_EP_3subBBCWest;
- TProfile2D *p_v2_VtxZsys[NEtaGaps];
- TProfile2D *p_v2_DCAsys[NEtaGaps];
- TProfile2D *p_v2_Nhitsys[NEtaGaps];
- TProfile2D *p_v3_VtxZsys[NEtaGaps];
- TProfile2D *p_v3_DCAsys[NEtaGaps];
- TProfile2D *p_v3_Nhitsys[NEtaGaps];*/
- // <vn> vs RunId & centrality - for RunQA
- // In this case, standard EP from 2-sub(TPC) method is used
- // TProfile2D *p_v2_vs_run;
- // TProfile2D *p_v3_vs_run;
- // p_v2cent_EP_3subTPC = new TProfile2D(Form("p_v2cent_EP_3subTPC"),Form("p_v2cent_EP_3subTPC"),16,-0.5,15.5,50,0.15,5.15);
- // p_v2cent_EP_3subBBCEast = new TProfile2D(Form("p_v2cent_EP_3subBBCEast"),Form("p_v2cent_EP_3subBBCEast"),16,-0.5,15.5,50,0.15,5.15);
- // p_v2cent_EP_3subBBCWest = new TProfile2D(Form("p_v2cent_EP_3subBBCWest"),Form("p_v2cent_EP_3subBBCWest"),16,-0.5,15.5,50,0.15,5.15);
- // p_v2_vs_run = new TProfile2D(Form("p_v2_vs_run"),Form("p_v2_vs_run"),fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- // p_v3_vs_run = new TProfile2D(Form("p_v3_vs_run"),Form("p_v3_vs_run"),fNBins, runId_min - 0.5, runId_max + 0.5, 16, -0.5, 15.5);
- for (int iGap = 0; iGap < NEtaGaps; iGap++)
- {
- for (int iCent = 0; iCent < 16; iCent++)
- {
- p_v2run_EP[iGap][iCent] = new TProfile2D(Form("p_v2_run_EP_gap%i_cent%i", iGap, iCent), Form("p_v2_run_EP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- p_v3run_EP[iGap][iCent] = new TProfile2D(Form("p_v3_run_EP_gap%i_cent%i", iGap, iCent), Form("p_v3_run_EP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- p_v2run_SP[iGap][iCent] = new TProfile2D(Form("p_v2_run_SP_gap%i_cent%i", iGap, iCent), Form("p_v2_run_SP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- p_v3run_SP[iGap][iCent] = new TProfile2D(Form("p_v3_run_SP_gap%i_cent%i", iGap, iCent), Form("p_v3_run_SP_gap%i_cent%i", iGap, iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- if (iGap == 0)
- {
- p_v2run_BBC[iCent] = new TProfile2D(Form("p_v2_run_BBC_cent%i", iCent), Form("p_v2_run_BBC_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- p_v2run_ZDC[iCent] = new TProfile2D(Form("p_v2_run_ZDC_cent%i", iCent), Form("p_v2_run_ZDC_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- p_v2run_ZDC_East[iCent] = new TProfile2D(Form("p_v2_run_ZDC_East_cent%i", iCent), Form("p_v2_run_ZDC_East_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- p_v2run_ZDC_West[iCent] = new TProfile2D(Form("p_v2_run_ZDC_West_cent%i", iCent), Form("p_v2_run_ZDC_West_cent%i", iCent), fNBins, runId_min - 0.5, runId_max + 0.5, 50, 0.15, 5.15);
- }
- }
- }
- for (int iGap = 0; iGap < NEtaGaps; iGap++)
- {
- p_v2cent_EP_ptgaps[iGap] = new TProfile2D(Form("p_v2cent_EP_gap%i_ptgaps", iGap), Form("p_v2cent_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- p_v3cent_EP_ptgaps[iGap] = new TProfile2D(Form("p_v3cent_EP_gap%i_ptgaps", iGap), Form("p_v3cent_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- p_v2cent_SP_ptgaps[iGap] = new TProfile2D(Form("p_v2cent_SP_gap%i_ptgaps", iGap), Form("p_v2cent_SP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- p_v3cent_SP_ptgaps[iGap] = new TProfile2D(Form("p_v3cent_SP_gap%i_ptgaps", iGap), Form("p_v3cent_SP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- p_v2WEST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v2WEST_EP_gap%i_ptgaps", iGap), Form("p_v2WEST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- p_v2EAST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v2EAST_EP_gap%i_ptgaps", iGap), Form("p_v2EAST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- p_v3WEST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v3WEST_EP_gap%i_ptgaps", iGap), Form("p_v3WEST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- p_v3EAST_EP_ptgaps[iGap] = new TProfile2D(Form("p_v3EAST_EP_gap%i_ptgaps", iGap), Form("p_v3EAST_EP_gap%i_ptgaps", iGap), 16, -0.5, 15.5, 50, 0.15, 5.15);
- /*p_v2_VtxZsys[iGap] = new TProfile2D(Form("p_v2_VtxZsys%i",iGap),Form("p_v2_VtxZsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
- p_v2_DCAsys[iGap] = new TProfile2D(Form("p_v2_DCAsys%i",iGap),Form("p_v2_DCAsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
- p_v2_Nhitsys[iGap] = new TProfile2D(Form("p_v2_Nhitsys%i",iGap),Form("p_v2_Nhitsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
- p_v3_VtxZsys[iGap] = new TProfile2D(Form("p_v3_VtxZsys%i",iGap),Form("p_v3_VtxZsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
- p_v3_DCAsys[iGap] = new TProfile2D(Form( "p_v3_DCAsys%i",iGap),Form( "p_v3_DCAsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);
- p_v3_Nhitsys[iGap] = new TProfile2D(Form("p_v3_Nhitsys%i",iGap),Form("p_v3_Nhitsys%i",iGap),16,-0.5,15.5,50,0.15,5.15);*/
- }
- for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- {
- p_vnBBC[iHarm] = new TProfile2D(Form("p_vnBBC%i", iHarm), Form("p_vnBBC%i", iHarm), 16, -0.5, 15.5, 50, 0.15, 5.15);
- // p_vnBBCEast[iHarm] = new TProfile2D(Form("p_vnBBCEast%i",iHarm),Form("p_vnBBCEast%i",iHarm),16,-0.5,15.5,50,0.15,5.15);
- // p_vnBBCWest[iHarm] = new TProfile2D(Form("p_vnBBCWest%i",iHarm),Form("p_vnBBCWest%i",iHarm),16,-0.5,15.5,50,0.15,5.15);
- }
- //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 fQv2FullEPCorr;
- TVector2 fQv2FullSP;
- TVector2 fQv3FullEP;
- TVector2 fQv3FullSP;
- TVector2 fQv2WestEP[NEtaGaps];
- TVector2 fQv2EastEP[NEtaGaps];
- TVector2 fQv2WestSP[NEtaGaps];
- TVector2 fQv2EastSP[NEtaGaps];
- TVector2 fQv3WestEP[NEtaGaps];
- TVector2 fQv3EastEP[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, dPsi;
- 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;
- Int_t shiftBin;
- Double_t psiSP_rec;
- Double_t psiEP_rec, psiEP_BBC_Full[fNharmonics + 1], shiftCos, shiftSin;
- Double_t psiEP_BBC_East[fNharmonics + 1], psiEP_BBC_West[fNharmonics + 1];
- // ZDC related EP determination
- Double_t psiEP_ZDC_Full[fNharmonics + 1], psiEP_ZDC_East[fNharmonics + 1], psiEP_ZDC_West[fNharmonics + 1];
- TVector2 qRecEastZDC[fNharmonics + 1], qRecWestZDC[fNharmonics + 1], qRecFullZDC[fNharmonics + 1];
- Double_t qxEast, qyEast, adcxEast, adcyEast;
- Double_t qxWest, qyWest, adcxWest, adcyWest;
- Double_t qxFull, qyFull, adcxFull, adcyFull;
- Double_t qxRecEast, qyRecEast;
- Double_t qxRecWest, qyRecWest;
- Double_t qxRecFull, qyRecFull;
- Double_t resEastWestZDC[fNharmonics + 1];
- // 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;
- }
- // Event selection
- if (!isGoodEvent(event))
- continue;
- goodEventCounter++;
- TVector3 pVtx = event->primaryVertex();
- // Track analysis
- Int_t nTracks = dst->numberOfTracks();
- //Counters
- fUsedTracks = 0;
- fQcountFull = 0;
- fQcountFullEast = 0;
- fQcountFullWest = 0;
- psiShiftedTMP = 0.;
- PsiSin = 0.;
- PsiCos = 0.;
- fQv2FullEP.Set(0., 0.);
- fQv2FullSP.Set(0., 0.);
- fQv3FullEP.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.);
- fQv2WestSP[i].Set(0., 0.);
- fQv2EastSP[i].Set(0., 0.);
- fQv3WestEP[i].Set(0., 0.);
- fQv3EastEP[i].Set(0., 0.);
- fQv3WestSP[i].Set(0., 0.);
- fQv3EastSP[i].Set(0., 0.);
- }
- // Vertex sign
- VtxSign = GetVzBin(pVtx.Z());
- if (VtxSign == -1)
- continue;
- //BBC EP
- 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 psiFull[fNharmonics + 1];
- TVector2 qGainFull[fNharmonics + 1];
- TVector2 qRecFull[fNharmonics + 1];
- Double_t psiEast[fNharmonics + 1];
- TVector2 qGainEast[fNharmonics + 1];
- TVector2 qRecEast[fNharmonics + 1];
- Double_t psiWest[fNharmonics + 1];
- TVector2 qGainWest[fNharmonics + 1];
- TVector2 qRecWest[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;
- //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();
- if (EventPlaneGainWest < 0.0)
- {
- EventPlaneGainWest += 2. * TMath::Pi() / (iHarm + 1.);
- }
- }
- psiWest[iHarm] = EventPlaneGainWest;
- qGainWest[iHarm] = QWestGain;
- //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();
- if (EventPlaneGainEast < 0.0)
- {
- EventPlaneGainEast += 2. * TMath::Pi() / (iHarm + 1.);
- }
- }
- psiEast[iHarm] = EventPlaneGainEast;
- qGainEast[iHarm] = QEastGain;
- }
- // 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()]);
- qRecEast[iHarm] = qGainEast[iHarm] - recEastTMP;
- qRecWest[iHarm] = qGainWest[iHarm] - recWestTMP;
- qRecFull[iHarm] = qGainFull[iHarm] - recFullTMP;
- }
- 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;
- 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.);
- // 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;
- 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.);
- // 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_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.);
- }
- //EP from BBC is ready: psiEP_BBC_Full[iHarm], where iHarm = 0 for 1st harmonics, 1 for 2nd & 2 for 3rd.
- // ZDC EP
- if (event->zdcSumAdcEast() < 1.)
- continue;
- if (event->zdcSumAdcWest() < 1.)
- continue;
- qxEast = 0.;
- qyEast = 0.;
- adcxEast = 0.;
- adcyEast = 0.;
- qxWest = 0.;
- qyWest = 0.;
- adcxWest = 0.;
- adcyWest = 0.;
- qxFull = 0.;
- qyFull = 0.;
- adcxFull = 0.;
- adcyFull = 0.;
- // Vertical ZDC-SMD
- for (int iStrip = 0; iStrip < fNZdcSmdStripsVertical; iStrip++)
- {
- qxEast += event->zdcSmdEastVertical(iStrip) * GetZDCPosition(0, 0, iStrip);
- adcxEast += event->zdcSmdEastVertical(iStrip);
- qxWest += event->zdcSmdWestVertical(iStrip) * GetZDCPosition(1, 0, iStrip);
- adcxWest += event->zdcSmdWestVertical(iStrip);
- }
- // Horizontal ZDC-SMD
- for (int iStrip = 0; iStrip < fNZdcSmdStripsHorizontal; iStrip++)
- {
- qyEast += event->zdcSmdEastHorizontal(iStrip) * GetZDCPosition(0, 1, iStrip);
- adcyEast += event->zdcSmdEastHorizontal(iStrip);
- qyWest += event->zdcSmdWestHorizontal(iStrip) * GetZDCPosition(1, 1, iStrip);
- adcyWest += event->zdcSmdWestHorizontal(iStrip);
- }
- if (adcxEast > 0)
- qxEast /= adcxEast;
- else
- qxEast = -999.;
- if (adcyEast > 0)
- qyEast /= adcyEast;
- else
- qyEast = -999.;
- if (adcxWest > 0)
- qxWest /= adcxWest;
- else
- qxWest = -999.;
- if (adcyWest > 0)
- qyWest /= adcyWest;
- else
- qyWest = -999.;
- // Full ZDC-SMD Q-vector
- if (qxEast == -999. || qxWest == -999.)
- qxFull = -999.;
- else
- qxFull = qxEast - qxWest;
- if (qyEast == -999. || qyWest == -999.)
- qyFull = -999.;
- else
- qyFull = qyEast - qyWest;
- if (qxEast != -999.)
- qxRecEast = qxEast - (Double_t)p_ZDC_Qx_East_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
- if (qyEast != -999.)
- qyRecEast = qyEast - (Double_t)p_ZDC_Qy_East_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
- if (qxWest != -999.)
- qxRecWest = qxWest - (Double_t)p_ZDC_Qx_West_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
- if (qyWest != -999.)
- qyRecWest = qyWest - (Double_t)p_ZDC_Qy_West_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
- if (qxFull != -999.)
- qxRecFull = qxFull - (Double_t)p_ZDC_Qx_Full_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
- if (qyFull != -999.)
- qyRecFull = qyFull - (Double_t)p_ZDC_Qy_Full_EP[0][VtxSign][(Int_t)event->cent16()][(Double_t)event->runId()];
- if (qxEast == -999.)
- qxRecEast = qxEast;
- if (qyEast == -999.)
- qyRecEast = qyEast;
- if (qxWest == -999.)
- qxRecWest = qxWest;
- if (qyWest == -999.)
- qyRecWest = qyWest;
- if (qxFull == -999.)
- qxRecFull = qxFull;
- if (qyFull == -999.)
- qyRecFull = qyFull;
- PsiSin = 0.;
- PsiCos = 0.;
- psiShiftedTMP = 0.;
- for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- {
- qRecEastZDC[iHarm].Set(qxRecEast, qyRecEast);
- qRecWestZDC[iHarm].Set(qxRecWest, qyRecWest);
- qRecFullZDC[iHarm].Set(qxRecFull, qyRecFull);
- // Full ZDC
- if (qRecFullZDC[iHarm].Mod())
- {
- psiEP_rec = qRecFullZDC[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_ZDC_Cos_Full_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- shiftSin = p_ZDC_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_ZDC_Full[iHarm] = psiEP_rec + dPsi;
- if (psiEP_ZDC_Full[iHarm] < 0.0)
- psiEP_ZDC_Full[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
- if (psiEP_ZDC_Full[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
- psiEP_ZDC_Full[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
- // East ZDC
- if (qRecEastZDC[iHarm].Mod())
- {
- psiEP_rec = qRecEastZDC[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_ZDC_Cos_East_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- shiftSin = p_ZDC_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_ZDC_East[iHarm] = psiEP_rec + dPsi;
- if (psiEP_ZDC_East[iHarm] < 0.0)
- psiEP_ZDC_East[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
- if (psiEP_ZDC_East[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
- psiEP_ZDC_East[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
- // West ZDC
- if (qRecWestZDC[iHarm].Mod())
- {
- psiEP_rec = qRecWestZDC[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_ZDC_Cos_West_EP[iHarm][VtxSign][iK][(Int_t)event->cent16()][(Double_t)event->runId()];
- shiftSin = p_ZDC_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_ZDC_West[iHarm] = psiEP_rec + dPsi;
- if (psiEP_ZDC_West[iHarm] < 0.0)
- psiEP_ZDC_West[iHarm] += 2. * TMath::Pi() / (iHarm + 1.);
- if (psiEP_ZDC_West[iHarm] > 2. * TMath::Pi() / (iHarm + 1.))
- psiEP_ZDC_West[iHarm] -= 2. * TMath::Pi() / (iHarm + 1.);
- }
- // Track loop
- for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
- {
- // Retrieve i-th femto track
- StFemtoTrack *femtoTrack = dst->track(iTrk);
- if (!femtoTrack)
- continue;
- // Must be a primary track
- if (!femtoTrack->isPrimary())
- continue;
- //Track selection
- if (!isGoodTrack(femtoTrack, energy, pVtx))
- continue;
- 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]++;
- }
- }
- // Check if track has TOF signal
- if (femtoTrack->isTofTrack())
- {
- } //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 = AngleShift(Psi_ReCenter + dPsi, 2.);
- }
- 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]));
- }
- }
- //FLOW CALCULATIONS
- // Second Track loop
- for (Int_t iTrk = 0; iTrk < nTracks; iTrk++)
- {
- // Retrieve i-th femto track
- StFemtoTrack *femtoTrack = dst->track(iTrk);
- if (!femtoTrack)
- continue;
- // Must be a primary track
- if (!femtoTrack->isPrimary())
- continue;
- //Track selection
- if (isGoodTrackFlow(femtoTrack, energy, pVtx))
- {
- weight = GetWeight(femtoTrack);
- // Flow with BBC EP
- for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- {
- res2EP = hResBBCFULL[iHarm]->GetBinContent(hResBBCFULL[iHarm]->FindBin(event->cent16()));
- flow2EP = TMath::Cos((iHarm + 1.) * (femtoTrack->phi() - psiEP_BBC_Full[iHarm])) / res2EP;
- p_vnBBC[iHarm]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
- if (iHarm == 1)
- {
- p_v2run_BBC[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
- }
- /*res2EP = hResBBCHalf[iHarm]->GetBinContent(hResBBCHalf[iHarm]->FindBin(event->cent16()));
- flow2EP = TMath::Cos((iHarm+1.)*(femtoTrack->phi()-psiEP_BBC_East[iHarm])) / res2EP;
- p_vnBBCEast[iHarm]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- res2EP = hResBBCHalf[iHarm]->GetBinContent(hResBBCHalf[iHarm]->FindBin(event->cent16()));
- flow2EP = TMath::Cos((iHarm+1.)*(femtoTrack->phi()-psiEP_BBC_West[iHarm])) / res2EP;
- p_vnBBCWest[iHarm]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);*/
- }
- // Flow with ZDC EP
- // Full
- res2EP = hResZDCFULL[1]->GetBinContent(hResZDCFULL[1]->FindBin(event->cent16()));
- flow2EP = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_Full[0])) / res2EP;
- p_v2run_ZDC[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
- // East
- res2EP = hResZDCHalf[1]->GetBinContent(hResZDCHalf[1]->FindBin(event->cent16()));
- flow2EP = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_East[0])) / res2EP;
- p_v2run_ZDC_East[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
- // West
- res2EP = hResZDCHalf[1]->GetBinContent(hResZDCHalf[1]->FindBin(event->cent16()));
- flow2EP = TMath::Cos(2. * (femtoTrack->phi() - psiEP_ZDC_West[0])) / res2EP;
- p_v2run_ZDC_West[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
- ///////////////////////////////////////////////////////////////////////////////////////////
- //p_ptrun[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), weight);
- //p1_ptrun[(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), 1.);
- /*
- // 3-sub method
-
- //Get rid of autocorrelations
- //Since pt<2 GeV/c cut was used in EP determination, only tracks within this range should be corrected
- if (isGoodTrack(femtoTrack, energy, pVtx))
- {
- QvTMP.SetX(p_Qx_FULL_EP[0][VtxSign]->GetBinContent(p_Qx_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent16())));
- QvTMP.SetY(p_Qy_FULL_EP[0][VtxSign]->GetBinContent(p_Qy_FULL_EP[0][VtxSign]->FindBin((Double_t)event->runId(),(Double_t)event->cent16())));
- fQv2FullEPCorr = fQv2FullEP - weight*(CalcQvector(femtoTrack,2) - QvTMP);
- QvTMP.Set(0.,0.);
- Psi_ReCenter = TMath::ATan2(fQv2FullEPCorr.Y(), fQv2FullEPCorr.X()) / 2.;
- dPsi = 0.;
- for (int iK = 0; iK < NShiftOrderMax; iK++)
- {
- bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t)event->runId(), (Double_t)event->cent16());
- mCos = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);
- bin = p_Cos2_FULL_EP[VtxSign][iK]->FindBin((Double_t)event->runId(), (Double_t)event->cent16());
- mSin = p_Cos2_FULL_EP[VtxSign][iK]->GetBinContent(bin);
- dPsi += (1. / 2.) * (2. / (iK + 1)) * (-1. * mSin * TMath::Cos(shiftOrderCoeff2.at(iK) * Psi_ReCenter) + mCos * TMath::Sin(shiftOrderCoeff2.at(iK) * Psi_ReCenter));
- }
- Psi2FullEP = AngleShift(Psi_ReCenter + dPsi, 2.);
- }
- res2EPsqr = p_res2_EP_3sub_AB->GetBinContent(p_res2_EP_3sub_AB->FindBin(event->cent16())) *
- p_res2_EP_3sub_AC->GetBinContent(p_res2_EP_3sub_AC->FindBin(event->cent16())) /
- p_res2_EP_3sub_BC->GetBinContent(p_res2_EP_3sub_BC->FindBin(event->cent16()));
- res2EP = TMath::Sqrt(res2EPsqr);
- flow2EP = TMath::Cos(2.*(femtoTrack->phi()-Psi2FullEP)) / res2EP;
- p_v2cent_EP_3subTPC->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- // Calculate flow
- res2EPsqr = p_res2_EP_3sub_AB->GetBinContent(p_res2_EP_3sub_AB->FindBin(event->cent16())) *
- p_res2_EP_3sub_BC->GetBinContent(p_res2_EP_3sub_BC->FindBin(event->cent16())) /
- p_res2_EP_3sub_AC->GetBinContent(p_res2_EP_3sub_AC->FindBin(event->cent16()));
- res2EP = TMath::Sqrt(res2EPsqr);
- flow2EP = TMath::Cos(2.*(femtoTrack->phi()-psiEP_BBC_East[1])) / res2EP;
- p_v2cent_EP_3subBBCEast->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- res2EPsqr = p_res2_EP_3sub_AC->GetBinContent(p_res2_EP_3sub_AC->FindBin(event->cent16())) *
- p_res2_EP_3sub_BC->GetBinContent(p_res2_EP_3sub_BC->FindBin(event->cent16())) /
- p_res2_EP_3sub_AB->GetBinContent(p_res2_EP_3sub_AB->FindBin(event->cent16()));
- res2EP = TMath::Sqrt(res2EPsqr);
- flow2EP = TMath::Cos(2.*(femtoTrack->phi()-psiEP_BBC_West[1])) / res2EP;
- p_v2cent_EP_3subBBCWest->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- ///////////////////////////////////////////////////////////////////////////////////////////
- */
- for (int iGap = 0; iGap < NEtaGaps; iGap++)
- {
- res2EPsqr = p_res2_EP[iGap]->GetBinContent(p_res2_EP[iGap]->FindBin(event->cent16()));
- res2EP = TMath::Sqrt(res2EPsqr);
- res3EPsqr = p_res3_EP[iGap]->GetBinContent(p_res3_EP[iGap]->FindBin(event->cent16()));
- res3EP = TMath::Sqrt(res3EPsqr);
- res2SPsqr = p_res2_SP[iGap]->GetBinContent(p_res2_SP[iGap]->FindBin(event->cent16()));
- if (res2SPsqr < 0)
- std::cout << "Event " << iEvent << " | Res2SP NEGATIVEE: " << res2SPsqr << std::endl;
- if (res2SPsqr == 0)
- std::cout << "Res2SP ZEROOOOO" << std::endl;
- res2SP = TMath::Sqrt(res2SPsqr);
- res3SPsqr = p_res3_SP[iGap]->GetBinContent(p_res3_SP[iGap]->FindBin(event->cent16()));
- res3SP = TMath::Sqrt(res3SPsqr);
- //EAST
- //if (track->eta() > -1.*cutEta && track->eta() < 0.)
- if (femtoTrack->eta() > -1. * cutEta && femtoTrack->eta() < -1. * cutEtaGap.at(iGap))
- {
- flow2EP = TMath::Cos(2. * (femtoTrack->phi() - Psi2WestEP[iGap])) / res2EP;
- flow3EP = TMath::Cos(3. * (femtoTrack->phi() - Psi3WestEP[iGap])) / res3EP;
- p_v2cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
- p_v3cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);
- p_v2EAST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
- p_v3EAST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);
- p_v2run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
- p_v3run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP, weight);
- /*if (pVtx.Z() < 20.)
- {
- p_v2_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- p_v3_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
- }
- if (femtoTrack->gDCA(pVtx).Mag() < 2.)
- {
- p_v2_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- p_v3_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
- }
- if (femtoTrack->nHitsFit() > 10.)
- {
- p_v2_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- p_v3_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
- }
- if (iGap == 0)
- {
- p_v2_vs_run->Fill(event->runId(),event->cent16(),flow2EP);
- p_v3_vs_run->Fill(event->runId(),event->cent16(),flow3EP);
- }*/
- fUvSP.Set(1. * TMath::Cos(2. * femtoTrack->phi()), 1. * TMath::Sin(2. * femtoTrack->phi()));
- flow2SP = fUvSP * fQv2WestSP_Shift[iGap] / res2SP;
- fUvSP.Set(1. * TMath::Cos(3. * femtoTrack->phi()), 1. * TMath::Sin(3. * femtoTrack->phi()));
- flow3SP = fUvSP * fQv3WestSP_Shift[iGap] / res3SP;
- p_v2cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2SP, weight);
- p_v3cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3SP, weight);
- p_v2run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2SP, weight);
- p_v3run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3SP, weight);
- }
- //WEST
- //if (femtoTrack->eta() < 1.*cutEta && femtoTrack->eta() > 0.)
- if (femtoTrack->eta() < 1. * cutEta && femtoTrack->eta() > 1. * cutEtaGap.at(iGap))
- {
- flow2EP = TMath::Cos(2. * (femtoTrack->phi() - Psi2EastEP[iGap])) / res2EP;
- flow3EP = TMath::Cos(3. * (femtoTrack->phi() - Psi3EastEP[iGap])) / res3EP;
- p_v2cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
- p_v3cent_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);
- p_v2WEST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2EP, weight);
- p_v3WEST_EP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3EP, weight);
- p_v2run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2EP, weight);
- p_v3run_EP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3EP, weight);
- /*if (pVtx.Z() < 20.)
- {
- p_v2_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- p_v3_VtxZsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
- }
- if (femtoTrack->gDCA(pVtx).Mag() < 2.)
- {
- p_v2_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- p_v3_DCAsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
- }
- if (femtoTrack->nHitsFit() > 10.)
- {
- p_v2_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow2EP,weight);
- p_v3_Nhitsys[iGap]->Fill(event->cent16(),femtoTrack->pt(),flow3EP,weight);
- }
- if (iGap == 0)
- {
- p_v2_vs_run->Fill(event->runId(),event->cent16(),flow2EP);
- p_v3_vs_run->Fill(event->runId(),event->cent16(),flow3EP);
- }*/
- fUvSP.Set(1. * TMath::Cos(2. * femtoTrack->phi()), 1. * TMath::Sin(2. * femtoTrack->phi()));
- flow2SP = fUvSP * fQv2EastSP_Shift[iGap] / res2SP;
- fUvSP.Set(1. * TMath::Cos(3. * femtoTrack->phi()), 1. * TMath::Sin(3. * femtoTrack->phi()));
- flow3SP = fUvSP * fQv3EastSP_Shift[iGap] / res3SP;
- p_v2cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow2SP, weight);
- p_v3cent_SP_ptgaps[iGap]->Fill(event->cent16(), femtoTrack->pt(), flow3SP, weight);
- p_v2run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow2SP, weight);
- p_v3run_SP[iGap][(Int_t)event->cent16()]->Fill(event->runId(), femtoTrack->pt(), flow3SP, weight);
- }
- }
- }
- }
- } //for(Long64_t iEvent=0; iEvent<events2read; iEvent++)
- femtoReader->Finish();
- TFile *output = new TFile(outFile, "recreate");
- for (int iGap = 0; iGap < NEtaGaps; iGap++)
- {
- // p_v2cent_EP_ptgaps[iGap]->Write();
- // p_v3cent_EP_ptgaps[iGap]->Write();
- // p_v2cent_SP_ptgaps[iGap]->Write();
- // p_v3cent_SP_ptgaps[iGap]->Write();
- p_v2EAST_EP_ptgaps[iGap]->Write();
- p_v3EAST_EP_ptgaps[iGap]->Write();
- p_v2WEST_EP_ptgaps[iGap]->Write();
- p_v3WEST_EP_ptgaps[iGap]->Write();
- for (int iCent = 0; iCent < 16; iCent++)
- {
- p_v2run_EP[iGap][iCent]->Write();
- p_v3run_EP[iGap][iCent]->Write();
- p_v2run_SP[iGap][iCent]->Write();
- p_v3run_SP[iGap][iCent]->Write();
- if (iGap == 0)
- {
- p_v2run_BBC[iCent]->Write();
- }
- }
- for (int iCent = 0; iCent < 16; iCent++)
- {
- if (iGap == 0)
- {
- p_v2run_ZDC[iCent]->Write();
- p_v2run_ZDC_East[iCent]->Write();
- p_v2run_ZDC_West[iCent]->Write();
- }
- }
- for (int iCent = 0; iCent < 16; iCent++)
- {
- if (iGap == 0)
- {
- //p_ptrun[iCent]->Write();
- //p1_ptrun[iCent]->Write();
- }
- }
- /*p_v2_VtxZsys[iGap]->Write();
- p_v2_DCAsys[iGap]->Write();
- p_v2_Nhitsys[iGap]->Write();
- p_v3_VtxZsys[iGap]->Write();
- p_v3_DCAsys[iGap]->Write();
- p_v3_Nhitsys[iGap]->Write();*/
- }
- for (int iHarm = 0; iHarm < fNharmonics + 1; iHarm++)
- {
- p_vnBBC[iHarm]->Write();
- //p_vnBBCEast[iHarm]->Write();
- //p_vnBBCWest[iHarm]->Write();
- }
- /*p_v2cent_EP_3subTPC->Write();
- p_v2cent_EP_3subBBCEast->Write();
- p_v2cent_EP_3subBBCWest->Write();
- p_v2_vs_run->Write();
- p_v3_vs_run->Write();*/
- output->Close();
- std::cout << "I'm done with analysis. We'll have a Nobel Prize, Master!"
- << std::endl;
- std::cout << "Good events: " << goodEventCounter << 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;
- }
- Bool_t isBadRun(StFemtoEvent *const &event)
- {
- Int_t run = (Int_t)event->runId();
- if (run == 12132026 ||
- run == 12134014 ||
- run == 12134051 ||
- run == 12135008 ||
- run == 12136006 ||
- run == 12136025 ||
- run == 12136047 ||
- run == 12136064)
- return true;
- else
- return false;
- }
- Double_t GetZDCPosition(Int_t eastwest, Int_t verthori, Int_t strip)
- // Get position of each slat;strip starts from 0
- {
- Double_t zdcsmd_x[7] = {0.5, 2, 3.5, 5, 6.5, 8, 9.5};
- Double_t zdcsmd_y[8] = {1.25, 3.25, 5.25, 7.25, 9.25, 11.25, 13.25, 15.25};
- Double_t mZDCSMDCenterex = 4.72466;
- Double_t mZDCSMDCenterey = 5.53629;
- Double_t mZDCSMDCenterwx = 4.39604;
- Double_t mZDCSMDCenterwy = 5.19968;
- if (eastwest == 0 && verthori == 0)
- return zdcsmd_x[strip] - mZDCSMDCenterex;
- if (eastwest == 1 && verthori == 0)
- return mZDCSMDCenterwx - zdcsmd_x[strip];
- if (eastwest == 0 && verthori == 1)
- return zdcsmd_y[strip] / sqrt(2.) - mZDCSMDCenterey;
- if (eastwest == 1 && verthori == 1)
- return zdcsmd_y[strip] / sqrt(2.) - mZDCSMDCenterwy;
- return -999.;
- }
|