#ifndef fiTQun_shared_h_seen #define fiTQun_shared_h_seen #include #include #include #include #include #include #include "TScatTableF.h" #include "hemi/table.h" #ifndef NOSKLIBRARIES #define nPMT_max 11146 #else #define nPMT_max 100000 #endif const static int nRngMX = 15; #define nPID 7 #define nTrkTyp 9 #define nexpbinm1 99999 const static int nEpsilonBins = 999; #define ntpdfmupar 7 #define ntpdfppar 10 #define npars1dcprof 9 #define nsect1dcprof 3 namespace fiTQun_parameters{ enum particleIndices {ig, ie, imu, ipip, ikp, ip, icg}; enum TrackIndices {itrkg, itrke, itrkmu, itrkpip, itrkkp, itrkp, itrkcg, itrkpi0, itrkpipknk}; } class fiTQun_shared { private: bool flgLoaded; // int fRngTyp[nPID]; bool fUseFitCprof; bool fUseFittpdf; TString fQconstdir; int SK_Ver; TString WCSim_PMTType; TString WCSim_Config; TString strSKVer; TVector3 TransDetToGlbl; TRotation RotDetToGlbl; double tpdfmurang[2]; double tpdfprang[nPID][2]; double tpdfpar[nPID][2][ntpdfmupar][ntpdfppar]; TGraphErrors *gtpar[nPID][2][ntpdfmupar]; int nR0bin,nth0bin,nmom[nPID],imomcur; double R0binwrecp,th0binwrecp,mombfrac; double *R0bins,*th0bins,*mombins[nPID]; float ***InTblcur[3];//cherenkov profile integral tables for current particle type float ***InTbl[nPID][3];//cherenkov profile integral tables TGraph *gNphot[nPID];//number of photons produced TGraph *gNphtmom[nPID]; TGraph *gSmax[nPID];//track length TGraph *gIniso[nPID][2]; TGraph *gPDGRange; //pdg range in water TGraph *gPDGTime; //time duration of track derived from pdgrange double Inisopar[nPID][2][(npars1dcprof+1)*nsect1dcprof+1]; double nphotpar[nPID][(npars1dcprof+1)*nsect1dcprof+1]; double smaxpar[nPID][(npars1dcprof+1)*nsect1dcprof+1]; double *InPar[nPID][3]; double **InParPtr[nPID]; double **InPtrCur; double momofst[nPID],momminstp[nPID],momminmax[nPID][2]; double lnmom[5]; int nsct2dRbin,nsct2dthbin; double sct2dRbinwrecp,sct2dthbinwrecp; double *sct2dRbins,*sct2dthbins; float *SctTbl2d[3]; int nsct3dRbinm1,nsct3dwlbinm1,nsct3dthbinm1; double sct3dRbinw,sct3dwlbinw,sct3dthbinw; double *sct3dRbins,*sct3dwlbins,*sct3dthbins; float ***SctTbl3d; static double expbinwrecp; static double expbins[nexpbinm1+1]; static double expTbl[nexpbinm1+1]; static double epsilon_binwrecp; static double epsilon_bins[nEpsilonBins+1]; static double epsilon_Tbl [nEpsilonBins+1]; double fconpoly(double xval, int npars, int nsect, double *par); void InterpolateIn(double *In, double R0, double costh0); void LoadProfiles(bool fFitCProf, bool fFittProf); void LoadCProf(int iPID); void Load1dpars(TFile *fin, const char *ctmp, int npartmp, double *par); void DeleteInstances(); void DeletePtclProf(int iPID); void DeleteProfiles(); fiTQun_shared( int aMaxNPmt = 11146 ); static int fMaxNPMT; static fiTQun_shared* staticthis; // Private variables previously stored in fiTQun class // Verbosity control int fQuiet; // Multi-segment muon fit int MSFitMethod; double MSElossMin; double MSThetRes; // Scattered light prediction int flgscat; // Time window definition: double twnd_def[2]; double PrefitVtxMinDwall; // Peak threshold double PkFlatThr; double thrTres[2]; bool flgMaskDeadTimePMT; int flgTVarCorr; double GainCorrFact; // Total charge constraint double ctotqpnl,cosmintotq; bool fmskbkPMT; int flgBiasCorr; bool flgtime; double nwtr; double darkrate; int flgdraw; int ofnmflg; char ofnmstr[256]; double chrgfact; double QEEff; static double Phi0[nPID]; static double dEdx[nPID]; double QmuConvFact; // Used to be totQFact double PIDCut_epip[2]; // Used to be PIDCut_mue double PIDCut_mupip[2]; int nPi0SeedScanPoints; double pi0PreFitVtxShiftThr; double pi0PreFitVtxShiftStepSize; int pi0PreFitVtxShiftNsteps; double RCCut_a[2]; double RCCut_1Ea0; double RCCut_1Mua0; //MS fit internal variables int nSegMax; double MSScattSig; bool flgfullscat; public: static fiTQun_shared* Get( int aMaxNPmt = 11146 ){ if (!staticthis) { staticthis = new fiTQun_shared( aMaxNPmt ); } return staticthis; } ~fiTQun_shared(); bool AreYouReady(){ return flgLoaded; } int MAXNPMT(){ return fMaxNPMT; } TScatTableF* top_scattable; TScatTableF* bot_scattable; TScatTableF* side_scattable; static const int PIDarr[nPID]; static const TString PIDnames[nPID]; static const double m0[nPID]; // static const double dEdx[nPID]; static const double pi0mass=134.9766; // +/- 0.0006 static const int nTrkPar[nTrkTyp]; static const int npar[3];//number of charge pdf parameters(Lo,Mid,Hi) static const int nSnglTrkParams = 8; static const int nPi0Params = 12; static const int nPipParams = 11; static const int nPDK_MuGammaParams = 14; double WAttenL[2];//Water attenuation length // Parameters for indirect light time PDF // 0: gamma // 1: sigma // 2: t0 static double inDirtPDFpars[3]; // LUT time PDFs bool flgTimePDFmode; void LoadLookupTableTimePDF(); static hemi::Table3D *tpdf_direct[nPID]; static hemi::Table3D *tpdf_indirect[nPID]; // Parameters for vertex pre-fit: static int npfitr; // Number of pre-fit iterations. Maximum is 10. static double sgmar[10]; // Width (2*sigma^2) of the gaussian for the vertex estimator static double grdszVtx; // Spatial step size for vertex grid search static double grdsztime; // Temporal step size for vertex grid search // Parameter for MR fit static int nRingMaxTrials; // Number of rings to fit with maximum number of ring hypothesis (further rings fit only with one ring hypothesis -- e-like) int PMTtype[nPMT_max]; double PMTQE[nPMT_max]; double PMTphi[nPMT_max]; double PMT_R[nPMT_max]; double PMTpos[nPMT_max][3]; double PMTdir[nPMT_max][3]; Float_t PMTRadius; double PMTmaxPos[3]; int DetGeomType;// Detector geometry type double DetGeomR;// Detector radius double DetGeomL;// Detector half length // Define Additional Hyper-K variables (indexed by quadrant) // coordinates taken from: http://indico.ipmu.jp/indico/getFile.py/access?contribId=45&sessionId=10&resId=0&materialId=slides&confId=10 // ID corner points (x,y) = (1203.8,2109.5), (-1203.8,2109.5), (1100.4,-2109.5), (-1100.4,-2109.5) // top circle center (x,y) and radius = (+/-800.0,0) and 2909.5 // bottom circle center (x,y) and radius = (+/-600.0,0) and 2709.5 double DetGeomHKR[4]; // Radius of curvature (by quadrant) double DetGeomHKArcXCenter[4]; // Center of radius of curvature (by quadrant) double DetGeomHKXMax; // Maximum x-value of ID wall double DetGeomHKYMax; // Maximum y-value of ID wall double DetGeomHKXMaxAtYMax; // Maximum x-value at maximum y-value (top of tank) double DetGeomHKXMaxAtYMin; // Maximum x-value at minimum y-value (bottom of tank) int writeHistos; //write debug histos to tree void SetWriteHistos(int w){ writeHistos = w;} int GetWriteHistos() {return writeHistos;} void ReadParameters(int fScat=1, int *fRngtmp=NULL, bool fFitCProf=true, bool fFittProf=true, int iDetGeom=1); TVector3 TfmDetToGlblCoord(TVector3 tmpVct, bool flgPos, bool flgInverse=false); int GetSKVer(){return SK_Ver;} TString GetWCSimConfig(){return WCSim_Config;} TString GetWCSimPMTType(){return WCSim_PMTType;} int ClrInfo(int iMod); void SetTypemom(int iPID, double mom, double *In, double& nphot, double& smax); void EvalIn(double *In, double R0, double costh0); void GetTimeCoeff(int iPID, double momtmp, double *tcmupar); void Boundlogmu(double& logmu); double pCherenkovThr(int iPID); double ConvNphot2mom(int iPID, double nphottmp){ if (GetRingType(iPID)) return gNphtmom[iPID]->Eval(nphottmp); else return 0.; } int GetRingType(int iPID){return fRngTyp[iPID];} static double explintrp(double x); static double epsilon(double coseta); void IntrpScat3Dto2D(double *dwall); double GetScatRatio2D(int istp, double R, double costh); // double GetScatRatio3D(double R, double dwall, double costh); void SaveTimeWindow(int itwnd, double *twindow); void SavePreFit(int itwnd, double* fitparams); void GetPreFit(int itwnd, double* fitparams, double *twindow=NULL); void ClearPeaks(int itwnd); void SavePeak(int itwnd, double peakt0, double peakiness); void GetPeaks(int itwnd, int& npeaks, double *peakt0, double *peakiness); void SaveSubEvtInfo(int itwnd, int ipeak, int nHPMT, double Totalq); void SaveDefaultSnglTrkFit(double* fitparams, int itwnd, int ipeak, int iPID, double totmu, double nll, int PCflg, int *n50=NULL, double *q50=NULL); void GetDefaultSnglTrkFit(double* fitparams, int itwnd, int ipeak, int iPID, double& totmu, double& nll, int& PCflg); void SaveTestSnglTrkFit(double* fitparams, int istage, int itwnd, int iPID, double totmu, double nll, int PCflg); void SaveDefaultPi0Fit(double* fitparams, int index, double totmu, double nll, int PCflg); void SaveTestPi0Fit(double* fitparams, int istage, double totmu, double nll, int PCflg); void ClearMRFit(); void SaveMRFit(int idFit, int nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], double totmu, double nll, int PCflg, int iFit=-1); int GetMRFit(int idFit, int &nringtmp, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], double &totmu, double &nll, int &PCflg, int iFit=-1); int MoveMRFit(int idFit, int iDest); void GetListMRFit(int &nMRFit, int *idList, double *nllList, int *iSort); //MSmuon fit void SaveMSFit(int idMRFit, int nseg, int* segIndx, int iFit, int *TypeArr, double ParArr[][fiTQun_shared::nSnglTrkParams], double totmu, double nll, int PCflg); double EvalPdgRange(double mom); double EvalPdgTime(double mom); void SaveDefaultPDK_MuGammaFit(double* fitparams, int index, double totmu, double nll, int PCflg); // Public methods previously in the fiTQun class void Cantyoubequiet(int tmp) {fQuiet=tmp;} int IsQuiet() {return fQuiet;} void SetMSFitMethod(int msfitflg) {MSFitMethod=msfitflg;} int GetMSFitMethod() {return MSFitMethod;} void SetMSElossMin(double elossmin){MSElossMin=elossmin;} double GetMSElossMin() {return MSElossMin;} void SetMSThetRes(double thetres) {MSThetRes=thetres;} double GetMSThetRes() {return MSThetRes;} void SetDefaultTimeWindow(double t_start, double t_end) {twnd_def[0]=t_start; twnd_def[1]=t_end;} double GetDefaultTimeWindow(int index) { return twnd_def[index]; } void SetScatflg(int flg) {flgscat=flg;} int GetScatflg(){ return flgscat;} void SetPrefitVtxMinDwall(double tmp) {PrefitVtxMinDwall=tmp;} double GetPrefitVtxMinDwall() {return PrefitVtxMinDwall;} void SetPeakThr(double ctmp) {PkFlatThr=ctmp;} double GetPeakThr() { return PkFlatThr;} void SetThrTres(double *Ttmp) {thrTres[0]=Ttmp[0]; thrTres[1]=Ttmp[1];} void SetThrTres(int index, double Ttmp) {thrTres[index]=Ttmp;} double GetThrTres(int index) { return thrTres[index]; } void SetMaskDeadTimePMT(bool tmp) {flgMaskDeadTimePMT=tmp;} bool GetMaskDeadTimePMT() {return flgMaskDeadTimePMT;} void SetTVarCorrflg(int tmp) {flgTVarCorr=tmp;} int GetTVarCorrflg() {return flgTVarCorr;} void SetGainCorrFact(double tmp) {GainCorrFact=tmp;} double GetGainCorrFact() {return GainCorrFact;} void SetTotqMaxAng(double ctmp);//maximum angle in degree // double GetTotqMaxAng(){return cosmintotq;} void SetCosMinTotQ(double ctmp){ cosmintotq = ctmp; } double GetCosMinTotQ() {return cosmintotq; } void SetFmskbkPMT( bool bool_tmp) { fmskbkPMT = bool_tmp; } bool GetFmskbkPMT() { return fmskbkPMT; } void SetBiasCorrflg(int tmp) {flgBiasCorr=tmp;} int GetBiasCorrflg() { return flgBiasCorr; } void SetTimeflg(bool flg) {flgtime=flg;} bool GetTimeflg(){ return flgtime; } void SetWAttL(double watltmp) {for (int i=0; i<2; i++) WAttenL[i]=watltmp;} void Setnwtr(double nwtrtmp) {nwtr=nwtrtmp;} double Getnwtr() { return nwtr; } void SetDarkRate(double ctmp) {darkrate=ctmp;} double GetDarkRate() { return darkrate; } void SetDrawflg(int flg) {flgdraw=flg;} int GetDrawflg() { return flgdraw; } void SetActivePMTs(int iflgMode=1); void Setofilenm(char *ofname); char* Getofilenm() {return ofnmstr;} void Setofileflg(int flg) {ofnmflg = flg;} int Getofileflg() {return ofnmflg;} double GetChrgFact() { return chrgfact; } void SetChrgFact(double fact) { chrgfact = fact; } void SetQEEff(double QEval) {QEEff=QEval;} double GetQEEff() { return QEEff; } void SetPhi0(int iPhi0, double Phi0val); double GetPhi0(int iPhi0){ if (iPhi0 >= 0 and iPhi0 < nPID) return Phi0[iPhi0]; else return -1.; } void SetdEdx(int iPID, double val); double GetdEdx(int iPhi0){ if (iPhi0 >= 0 and iPhi0 < nPID) return dEdx[iPhi0]; else return -1.; } void SetQmuConvFact(double ctmp) {QmuConvFact=ctmp;} double GetQmuConvFact() {return QmuConvFact; } void SetPIDCutEPip(double tmp0=-20., double tmp1=0.) {PIDCut_epip[0]=tmp0; PIDCut_epip[1]=tmp1;} double GetPIDCutEPip(int index) {return PIDCut_epip[index];} void SetPIDCutMuPip(double tmp0=20., double tmp1=0.2) {PIDCut_mupip[0]=tmp0; PIDCut_mupip[1]=tmp1;} double GetPIDCutMuPip(int index) {return PIDCut_mupip[index];} void SetPi0PreFitVtxShiftThr(double thr){pi0PreFitVtxShiftThr = thr;} double GetPi0PreFitVtxShiftThr(){return pi0PreFitVtxShiftThr;} void SetPi0PreFitVtxShiftStepSize(double step){pi0PreFitVtxShiftStepSize = step;} double GetPi0PreFitVtxShiftStepSize(){return pi0PreFitVtxShiftStepSize;} void SetPi0PreFitVtxShiftNsteps(int n){pi0PreFitVtxShiftNsteps = n;} int GetPi0PreFitVtxShiftNsteps(){return pi0PreFitVtxShiftNsteps;} void SetNPi0SeedScanPoints(int npts){nPi0SeedScanPoints = npts;} int GetNPi0SeedScanPoints(){return nPi0SeedScanPoints;} void SetRCCut(double tmpa0=120., double tmpa1=-0.6, double tmp1Ea0=150., double tmp1Mua0=150.) { RCCut_a[0]=tmpa0; RCCut_a[1]=tmpa1; RCCut_1Ea0=tmp1Ea0; RCCut_1Mua0=tmp1Mua0; } double GetRCCut(int index){ return RCCut_a[index]; } double GetRCCutE() {return RCCut_1Ea0;} double GetRCCutMu() {return RCCut_1Mua0;} void SetNSegMax(int nsegmax) {nSegMax=nsegmax;} int GetNSegMax(){return nSegMax;} void SetMSScattSig(double scattsig) {MSScattSig=scattsig;} double GetMSScattSig() {return MSScattSig;} void SetTotqCoef(double ctmp) {ctotqpnl=ctmp;} double GetTotqCoef() {return ctotqpnl;} void SetFullScatflg(bool btemp) {flgfullscat = btemp;} bool GetFullScatflg(){ return flgfullscat;} void SetTimePDFMode(bool flg); bool GetTimePDFMode() {return flgTimePDFmode;} // Variables used in likelihood loop that need to be shared between fitters (static for performance) static int flgHit[nPMT_max]; static int flgMask[nPMT_max]; static int ihitPMT[nPMT_max]; static int nhitPMT; static double tsigmvtx; static double chrg[nPMT_max]; static double mutot; static double mu[nRngMX][nPMT_max]; static double muscat[nRngMX][nPMT_max]; static double tHit[nPMT_max]; static double tau[nRngMX][nPMT_max]; static double nglnLarr[nPMT_max]; static double nglnLtarr[nPMT_max]; static double costheta0[nPMT_max]; static double qtot_fwd; static double mutot_fwd; static int fRngTyp[nPID]; static int nring; static int PIDring[nRngMX];//Particle type of each ring }; #endif