00001 #include "NCUtils/NCOscProb.h"
00002 #include "NCUtils/NCType.h"
00003
00004
00005 #include "OscProb/OscCalc.h"
00006
00007 #include "MessageService/MsgService.h"
00008 #include "TComplex.h"
00009 #include "TMath.h"
00010
00011 #include <cassert>
00012 #include <cmath>
00013
00014 CVSID("$Id: NCOscProb.cxx,v 1.44 2014/02/18 03:18:27 rhatcher Exp $");
00015
00016 using namespace NCType;
00017 using std::vector;
00018 using NC::Utility::SQR;
00019
00020
00021 bool NC::OscProb::OscPars::IsEquiv(const NC::OscProb::OscPars* rhs) const
00022 {
00023 if(!rhs) return false;
00024 if(rhs->OscillationModel() != OscillationModel()) return false;
00025 for(int n = 0; n < NCType::kNumParameters; ++n){
00026 if(rhs->fVals[n].Uninitialized() != fVals[n].Uninitialized())
00027 return false;
00028 if(!fVals[n].Uninitialized() && rhs->fVals[n] != fVals[n])
00029 return false;
00030 }
00031 return true;
00032 }
00033
00034
00035 std::ostream& NC::OscProb::OscPars::Print(std::ostream& os) const
00036 {
00037 os << "[ ";
00038 for(int n = 0; n < NCType::kNumParameters; ++n)
00039 if(!fVals[n].Uninitialized()) os << fVals[n] << " ";
00040 os << "]";
00041 return os;
00042 }
00043
00044
00045 bool NC::OscProb::OscPars::operator<(const NC::OscProb::OscPars& rhs) const
00046 {
00047 if(IsEquiv(&rhs)) return false;
00048
00049 if(fModel < rhs.fModel) return true;
00050 if(fModel > rhs.fModel) return false;
00051
00052 for(int n = 0; n < NCType::kNumParameters; ++n){
00053 if(fVals[n].Uninitialized()) continue;
00054 assert(!rhs.fVals[n].Uninitialized());
00055 if(fVals[n] < rhs.fVals[n]) return true;
00056 if(fVals[n] > rhs.fVals[n]) return false;
00057 }
00058
00059 assert(0 && "Not reached");
00060 }
00061
00062
00063 std::ostream& NC::OscProb::operator<<(std::ostream& os,
00064 const NC::OscProb::OscPars& rhs)
00065 {
00066 return rhs.Print(os);
00067 }
00068
00069
00070 double NC::OscProb::OscProbFs(double ssDeltaMSqr31,
00071 double amu,
00072 double fs,
00073 double ae,
00074 int oscMode)
00075 {
00076 switch(oscMode){
00077 case NCType::kNuMuToNuMu: return 1. - amu*ssDeltaMSqr31;
00078 case NCType::kNuMuToNuS: return fs*amu*ssDeltaMSqr31;
00079 case NCType::kNuMuToNuE: return ae*ssDeltaMSqr31;
00080 case NCType::kNuMuToNuTau: return (amu*(1.-fs)-ae)*ssDeltaMSqr31;
00081 case NCType::kNuEToNuE: return 1;
00082 default:
00083 assert(0 && "Bad oscMode");
00084 }
00085 }
00086
00087
00088 double NC::OscProb::ThreeFlavor::
00089 TransitionProbability(NCType::EOscMode oscMode,
00090 NCType::EEventType interaction,
00091 double baseline,
00092 double energy) const
00093 {
00094 if(interaction == NCType::kNC){
00095 switch(oscMode){
00096 case NCType::kNuMuToNuTau: return 0.;
00097 case NCType::kNuMuToNuE: return 0.;
00098 case NCType::kNuMuToNuMu: return 1.;
00099 case NCType::kNuMuToNuS:
00100 assert(0 && "TODO - Old code would return 1 here, that's bad");
00101 case NCType::kNuEToNuE:
00102 return 1;
00103 }
00104 }
00105
00106 double pars[9];
00107 pars[OscPar::kTh12] = Theta12();
00108 pars[OscPar::kTh23] = Theta23();
00109 pars[OscPar::kTh13] = Theta13();
00110 pars[OscPar::kDeltaM23] = DeltaMSqr32()*Hierarchy();
00111 pars[OscPar::kDeltaM12] = DeltaMSqr12();
00112 pars[OscPar::kDelta] = Delta13();
00113 pars[OscPar::kDensity] = RockDensity();
00114 pars[OscPar::kL] = baseline;
00115 pars[OscPar::kNuAntiNu] = -1;
00116
00117
00118
00119 static OscCalc oc;
00120 oc.SetOscParam(pars);
00121
00122 switch(oscMode){
00123 case NCType::kNuMuToNuMu: return oc.Oscillate(14, 14, energy);
00124 case NCType::kNuMuToNuTau: return oc.Oscillate(16, 14, energy);
00125 case NCType::kNuMuToNuE: return oc.Oscillate(12, 14, energy);
00126 case NCType::kNuEToNuE: return oc.Oscillate(12, 12, energy);
00127 case NCType::kNuMuToNuS:
00128 assert(0 && "TODO - Old code would return 1 here - bad");
00129 }
00130
00131 assert(0 && "Not reached");
00132 }
00133
00134
00135 double NC::OscProb::SterileFraction::TransitionProbability(
00136 NCType::EOscMode oscMode,
00137 NCType::EEventType interaction,
00138 double baseline,
00139 double energy) const
00140 {
00141 double amu = 4.*UMu3Sqr()*(1.-UMu3Sqr());
00142 double ssDeltaMSqr32 = FindSinSqrDeltaMSqr(energy, baseline, DeltaMSqr32());
00143 double fs = Fs();
00144 double ae = 4.*UMu3Sqr()*UE3Sqr();
00145
00146 if(interaction == NCType::kNC){
00147 switch(oscMode){
00148 case NCType::kNuMuToNuTau: return 0.;
00149 case NCType::kNuMuToNuE: return 0.;
00150 case NCType::kNuEToNuE: return 1.;
00151 case NCType::kNuMuToNuS:
00152 return OscProbFs(ssDeltaMSqr32, amu, fs, ae,
00153 NCType::kNuMuToNuS);
00154 case NCType::kNuMuToNuMu:
00155 return 1.-OscProbFs(ssDeltaMSqr32, amu, fs, ae,
00156 NCType::kNuMuToNuS);
00157 }
00158 }
00159
00160 return OscProbFs(ssDeltaMSqr32, amu, fs, ae, oscMode);
00161 }
00162
00163
00164 double NC::OscProb::SterileFractionTauNorm::TransitionProbability(
00165 NCType::EOscMode oscMode,
00166 NCType::EEventType interaction,
00167 double baseline,
00168 double energy) const
00169 {
00170 double amu = 4.*UMu3Sqr()*(1.-UMu3Sqr());
00171 double ssDeltaMSqr32 = FindSinSqrDeltaMSqr(energy, baseline, DeltaMSqr32());
00172 double fs = Fs();
00173 double ae = 4.*UMu3Sqr()*UE3Sqr();
00174
00175 if(interaction == NCType::kNC){
00176 switch(oscMode){
00177 case NCType::kNuMuToNuTau: return 0.;
00178 case NCType::kNuMuToNuE: return 0.;
00179 case NCType::kNuEToNuE: return 1.;
00180 case NCType::kNuMuToNuS:
00181 return OscProbFs(ssDeltaMSqr32, amu, fs, ae,
00182 NCType::kNuMuToNuS);
00183 case NCType::kNuMuToNuMu:
00184 return 1.-OscProbFs(ssDeltaMSqr32, amu, fs, ae,
00185 NCType::kNuMuToNuS);
00186 }
00187 }
00188
00189 double ret=OscProbFs(ssDeltaMSqr32, amu, fs, ae, oscMode);
00190 if (oscMode==NCType::kNuMuToNuTau) ret*=TauScale();
00191
00192 if (ret<0) return 0;
00193 return ret;
00194 }
00195
00196
00197 double NC::OscProb::FindSinSqrDeltaMSqr(double E, double L, double dmsq)
00198 {
00199 return SQR(TMath::Sin(NCType::k127*L/E*dmsq));
00200 }
00201
00202
00203 NCType::EOscMode NC::OscProb::ToOscMode(int from, int to)
00204 {
00205
00206 if(from < 0) from *= -1;
00207 if(to < 0) to *= -1;
00208
00209
00210 if(from%2) ++from;
00211 if(to%2) ++to;
00212
00213 if(from == 14 && to == 14) return NCType::kNuMuToNuMu;
00214 if(from == 14 && to == 16) return NCType::kNuMuToNuTau;
00215 if(from == 14 && to == 12) return NCType::kNuMuToNuE;
00216 if(from == 14 && to == 0) return NCType::kNuMuToNuS;
00217 if(from == 12 && to == 12) return NCType::kNuEToNuE;
00218
00219 assert(0 && "Don't have an enum for this oscillation mode");
00220 }
00221
00222
00223 double NC::OscProb::FourFlavorBase::
00224 TransitionProbability(NCType::EOscMode mode,
00225 NCType::EEventType intType,
00226 double baseline,
00227 double trueE) const
00228 {
00229 if(intType == NCType::kNC){
00230 switch(mode){
00231 case NCType::kNuMuToNuMu:
00232 return 1-TransitionProbability(NCType::kNuMuToNuS, intType, baseline, trueE);
00233 case NCType::kNuMuToNuTau:
00234 case NCType::kNuMuToNuE:
00235 return 0;
00236 case NCType::kNuEToNuE:
00237 return 1;
00238 case NCType::kNuMuToNuS:
00239
00240 break;
00241 default:
00242 assert(0 && "Unknown mode");
00243 }
00244 }
00245
00246 fLoverE = baseline/trueE;
00247
00248
00249 const double s13 = TMath::Sin(Theta13());
00250 const double s23 = TMath::Sin(Theta23());
00251 const double s14 = TMath::Sin(Theta14());
00252 const double s24 = TMath::Sin(Theta24());
00253 const double s34 = TMath::Sin(Theta34());
00254
00255 const double c13 = TMath::Cos(Theta13());
00256 const double c23 = TMath::Cos(Theta23());
00257 const double c14 = TMath::Cos(Theta14());
00258 const double c24 = TMath::Cos(Theta24());
00259 const double c34 = TMath::Cos(Theta34());
00260
00261
00262 const TComplex emd1(1, -Delta1(), true);
00263 const TComplex epd1 = TComplex::Conjugate(emd1);
00264 const TComplex emd2(1, -Delta2(), true);
00265 const TComplex epd2 = TComplex::Conjugate(emd2);
00266
00267
00268 const TComplex ue3 = emd1*c14*s13;
00269
00270 const TComplex ue4 = s14;
00271
00272 const TComplex umu3 = emd1*emd2*-s13*s14*s24 + c13*s23*c24;
00273
00274 const TComplex umu4 = emd2*c14*s24;
00275
00276 const TComplex utau3_1 = emd1*-s13*s14*c24*s34;
00277 const TComplex utau3_2 = epd2*-c13*s23*s24*s34;
00278 const TComplex utau3 = c13*c23*c34 + utau3_1 + utau3_2;
00279
00280 const TComplex utau4 = c14*c24*s34;
00281
00282 const TComplex us3_1 = emd1*-s13*s14*c24*c34;
00283 const TComplex us3_2 = epd2*-c13*s23*s24*c34;
00284 const TComplex us3 = -c13*c23*s34 + us3_1 + us3_2;
00285
00286 const TComplex us4 = c14*c24*c34;
00287
00288
00289 const double ue3Sqr = ue3.Rho2();
00290 const double ue4Sqr = ue4.Rho2();
00291 const double umu3Sqr = umu3.Rho2();
00292 const double umu4Sqr = umu4.Rho2();
00293 const double utau3Sqr = utau3.Rho2();
00294 const double utau4Sqr = utau4.Rho2();
00295 const double us3Sqr = us3.Rho2();
00296 const double us4Sqr = us4.Rho2();
00297
00298
00299 const double epsilon = 1e-5;
00300 const double col3Sum = ue3Sqr+umu3Sqr+utau3Sqr+us3Sqr;
00301 const double col4Sum = ue4Sqr+umu4Sqr+utau4Sqr+us4Sqr;
00302 if(TMath::Abs(col3Sum-1) > epsilon ||
00303 TMath::Abs(col4Sum-1) > epsilon)
00304 MSG("NCOscProb", Msg::kWarning) << "Warning: unitarity violated!!!!"
00305 << endl;
00306
00307 const TComplex umu4conj = TComplex::Conjugate(umu4);
00308
00309
00310 const TComplex umuue = umu4conj*ue4 *umu3*TComplex::Conjugate(ue3);
00311 const TComplex umuutau = umu4conj*utau4*umu3*TComplex::Conjugate(utau3);
00312 const TComplex umuus = umu4conj*us4 *umu3*TComplex::Conjugate(us3);
00313
00314 double p = 0;
00315
00316 #ifdef CASEMUTOALPHA
00317 #error Someone already defined CASEMUTOALPHA, have to choose a different name
00318 #endif
00319
00320
00321 #define CASEMUTOALPHA(alpha) \
00322 p = umu3Sqr*u##alpha##3Sqr*sinSqrDelta31(); \
00323 p += umu4Sqr*u##alpha##4Sqr*sinSqrDelta41(); \
00324 p += (umuu##alpha).Re()*(sinSqrDelta31() - sinSqrDelta43() + sinSqrDelta41()); \
00325 p *= 4; \
00326 p += 2.*(umuu##alpha).Im()*(sin2Delta31() - sin2Delta41() + sin2Delta43()); \
00327 return p
00328
00329 switch(mode){
00330 case NCType::kNuMuToNuMu:
00331 p = umu3Sqr*(1-umu3Sqr-umu4Sqr)*sinSqrDelta31();
00332 p += umu4Sqr*(1-umu3Sqr-umu4Sqr)*sinSqrDelta41();
00333 p += umu4Sqr*umu3Sqr*sinSqrDelta43();
00334 p *= 4.;
00335 return 1 - p;
00336
00337 case NCType::kNuMuToNuTau:
00338 CASEMUTOALPHA(tau);
00339
00340 case NCType::kNuMuToNuE:
00341 CASEMUTOALPHA(e);
00342
00343 case NCType::kNuMuToNuS:
00344 CASEMUTOALPHA(s);
00345
00346 case NCType::kNuEToNuE:
00347
00348 return 1;
00349
00350 default:
00351 assert(0 && "Unknown mode");
00352 }
00353
00354 #undef CASEMUTOALPHA
00355 }
00356
00357
00358 double NC::OscProb::Decay::TransitionProbability(EOscMode mode,
00359 EEventType intType,
00360 double L,
00361 double E) const
00362 {
00363 const double ssq = SQR(TMath::Sin(Theta()));
00364 const double csq = SQR(TMath::Cos(Theta()));
00365 const double arg = L/(2*E);
00366 const double exp = TMath::Exp(-Alpha()*arg);
00367
00368
00369 const double dmsq = 4*NCType::k127*DeltaMSqr();
00370
00371 switch(mode){
00372 case kNuMuToNuMu:
00373 if(intType == NCType::kNC) return 1-(1-SQR(exp))*ssq;
00374
00375 return (SQR(csq) + SQR(ssq*exp) +
00376 2*csq*ssq*exp*TMath::Cos(dmsq*arg));
00377
00378 case kNuMuToNuTau:
00379 if(intType == NCType::kNC) return 0;
00380
00381 return csq*ssq*(1+SQR(exp)-2*TMath::Cos(dmsq*arg)*exp);
00382
00383 case kNuMuToNuE:
00384 return 0;
00385
00386
00387 case kNuMuToNuS:
00388 return (1-SQR(exp))*ssq;
00389
00390 case kNuEToNuE:
00391 return 1;
00392 }
00393
00394 assert(0 && "Unknown oscillation mode");
00395 }
00396
00397
00398 double NC::OscProb::Decoherence::TransitionProbability(EOscMode mode,
00399 EEventType intType,
00400 double baseline,
00401 double trueEnergy) const
00402 {
00403 switch(mode){
00404 case kNuMuToNuMu:
00405 if(intType == NCType::kNC) return 1;
00406 return 1-SQR(TMath::Sin(2*Theta()))/2*(1-TMath::Exp(-SQR(Mu())*baseline/(2*trueEnergy)));
00407
00408 case kNuMuToNuTau:
00409 return 0;
00410
00411 case kNuMuToNuE:
00412 return 0;
00413
00414 case kNuMuToNuS:
00415 return 0;
00416
00417 case kNuEToNuE:
00418 return 1;
00419 }
00420
00421 assert(0 && "Unknown oscillation mode");
00422 }