loon: OscCalc Class Reference

OscCalc::OscCalc ( )

Definition at line 12 of file DataUtil/OscCalc.cxx.

00013 {
00014   for(int i = 0; i < 9; i++) fOscPar[i] = 0;  
00015   double par[9] = {0};
00016   SetOscParam(par, true);
00017 }

OscCalc::OscCalc ( )

double OscCalc::ElecToElec ( double * x )

Definition at line 546 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References MuELoss::e, fCos2Th, fOscPar, fSin2Th, fSinDCP, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, and OscPar::kTh13.

00547 {
00548   double E = x[0]; //energy
00549                                                                                                                        
00550   //Building the standard terms
00551   double L = fOscPar[OscPar::kL];  //baseline
00552   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00553   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00554   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00555                                                                                                                        
00556   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00557   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00558                                                                                                                        
00559   double sin_th12 = fSinTh[OscPar::kTh12];
00560  
00561 //  double cos_2th23 = fCos2Th[OscPar::kTh23];
00562   double cos_2th13 = fCos2Th[OscPar::kTh13];
00563   double cos_2th12 = fCos2Th[OscPar::kTh12];
00564                                                                                                                        
00565   double d_cp = fOscPar[OscPar::kDelta];
00566   double sin_dcp = fSinDCP;
00567                                                                                                                        
00568   //Building the more complicated terms
00569   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00570   double A = 2*fV*E*1e9/(dmsq_13);
00571   double alpha = dmsq_12/dmsq_13;
00572                                                                                                                        
00573   // A and d_cp both change sign for antineutrinos
00574   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00575   A *= plusminus;
00576   d_cp *= plusminus;
00577   sin_dcp *= plusminus;
00578                                                                                                                        
00579   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00580   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00581                                                                                                                        
00582   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00583   if(fabs(alpha) > 1e-10){  //want to be careful here
00584     double temp = cos_2th12 - A/alpha;
00585     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00586   }
00587                                                                                                                        
00588   //More complicated sin and cosine terms
00589   double cosC13Delta = cos(C13*Delta);
00590   double sinC13Delta = sin(C13*Delta);
00591                                                                                                                        
00592   double cosaC12Delta = 0;
00593   double sinaC12Delta = 0;
00594                                                                                                                        
00595   if(fabs(alpha) > 1e-10){
00596     cosaC12Delta = cos(alpha*C12*Delta);
00597     sinaC12Delta = sin(alpha*C12*Delta);
00598   } // otherwise not relevant
00599                                                                                                                        
00600   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00601   // this is the equivalent of Eq 45 & 46 corrected for Mu to E instead of E to Mu
00602                                                                                                                        
00603   // Leading order term
00604   double p1 = 1 - sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00605                                                                                                                        
00606   // terms that appear at order alpha
00607   double p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 -
00608                       A*sinC13Delta*(cos_2th13-A)/(C13*C13);
00609                                                                                                                        
00610   double p2 = +2*sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;
00611   //  p1 + p2 is the complete contribution for this expansion
00612                                                                                                                        
00613   // Now for the expansion in orders of sin(th13) (good to all order alpha)
00614   //  this is the equivalent of Eq 63 and 64
00615                                                                                                                        
00616   // leading order term
00617   double pa1 = 1.0, pa2 = 0.0;
00618                                                                                                                        
00619   if(fabs(alpha) > 1e-10){
00620     // leading order term
00621     pa1 = 1.0 - sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);
00622   }
00623   //pa1 is the complete contribution from this expansion, there is no order s13^1 term
00624                                                                                                                        
00625   // In order to combine the information correctly we need to add the two
00626   //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
00627   //  these may be taken from the expansion to second order in both parameters
00628   //  Equation 30
00629                                                                                                                        
00630   double repeated = 1;
00631                                                                                                                        
00632   //  Calculate the total probability
00633   double totalP = p1+p2 + (pa1+pa2) - repeated;
00634                                                                                                                        
00635   return totalP;
00636 }

double OscCalc::ElecToElec ( double x )

Definition at line 476 of file DataUtil/OscCalc.cxx.

References MuELoss::e, fCos2Th, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, and OscPar::kTh13.

Referenced by Oscillate(), OscillateMatter(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), and Extrapolate2D::SetNominalOscProb().

00477 {
00478   double E = x; //energy
00479                                                                                                                        
00480   //Building the standard terms
00481   double L = fOscPar[OscPar::kL];  //baseline
00482   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00483   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00484   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00485                                                                                                                        
00486   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00487   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00488                                                                                                                        
00489   double sin_th12 = fSinTh[OscPar::kTh12];
00490  
00491 //  double cos_2th23 = fCos2Th[OscPar::kTh23];
00492   double cos_2th13 = fCos2Th[OscPar::kTh13];
00493   double cos_2th12 = fCos2Th[OscPar::kTh12];
00494                                                                                                                        
00495   double d_cp = fOscPar[OscPar::kDelta];
00496   double sin_dcp = fSinDCP;
00497                                                                                                                        
00498   //Building the more complicated terms
00499   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00500   double A = 2*fV*E*1e9/(dmsq_13);
00501   double alpha = dmsq_12/dmsq_13;
00502                                                                                                                        
00503   // A and d_cp both change sign for antineutrinos
00504   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00505   A *= plusminus;
00506   d_cp *= plusminus;
00507   sin_dcp *= plusminus;
00508                                                                                                                        
00509   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00510   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00511                                                                                                                        
00512   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00513   if(fabs(alpha) > 1e-10){  //want to be careful here
00514     double temp = cos_2th12 - A/alpha;
00515     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00516   }
00517                                                                                                                        
00518   //More complicated sin and cosine terms
00519   double cosC13Delta = cos(C13*Delta);
00520   double sinC13Delta = sin(C13*Delta);
00521                                                                                                                        
00522   double cosaC12Delta = 0;
00523   double sinaC12Delta = 0;
00524                                                                                                                        
00525   if(fabs(alpha) > 1e-10){
00526     cosaC12Delta = cos(alpha*C12*Delta);
00527     sinaC12Delta = sin(alpha*C12*Delta);
00528   } // otherwise not relevant
00529                                                                                                                        
00530   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00531   // this is the equivalent of Eq 45 & 46 corrected for Mu to E instead of E to Mu
00532                                                                                                                        
00533   // Leading order term
00534   double p1 = 1 - sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00535                                                                                                                        
00536   // terms that appear at order alpha
00537   double p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 -
00538                       A*sinC13Delta*(cos_2th13-A)/(C13*C13);
00539                                                                                                                        
00540   double p2 = +2*sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;
00541   //  p1 + p2 is the complete contribution for this expansion
00542                                                                                                                        
00543   // Now for the expansion in orders of sin(th13) (good to all order alpha)
00544   //  this is the equivalent of Eq 63 and 64
00545                                                                                                                        
00546   // leading order term
00547   double pa1 = 1.0, pa2 = 0.0;
00548                                                                                                                        
00549   if(fabs(alpha) > 1e-10){
00550     // leading order term
00551     pa1 = 1.0 - sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);
00552   }
00553   //pa1 is the complete contribution from this expansion, there is no order s13^1 term
00554                                                                                                                        
00555   // In order to combine the information correctly we need to add the two
00556   //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
00557   //  these may be taken from the expansion to second order in both parameters
00558   //  Equation 30
00559                                                                                                                        
00560   double repeated = 1;
00561                                                                                                                        
00562   //  Calculate the total probability
00563   double totalP = p1+p2 + (pa1+pa2) - repeated;
00564                                                                                                                        
00565   return totalP;
00566 }

double OscCalc::ElecToMu ( double * x )

Definition at line 528 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References fOscPar, fSinDCP, OscPar::kDelta, and MuToElec().

00529 {
00530   // Flip delta to reverse direction
00531   double oldSinDelta = fSinDCP;
00532   double oldDelta =  fOscPar[OscPar::kDelta];
00533 
00534   fOscPar[OscPar::kDelta] = -oldDelta;
00535   fSinDCP = -oldSinDelta;
00536                                                                                                                        
00537   double prob = OscCalc::MuToElec(x);
00538                                                                                                                        
00539   //restore the world
00540   fOscPar[OscPar::kDelta] = oldDelta;
00541   fSinDCP = oldSinDelta;
00542  
00543   return prob;   
00544 }

double OscCalc::ElecToMu ( double x )

Definition at line 458 of file DataUtil/OscCalc.cxx.

References fOscPar, fSinDCP, OscPar::kDelta, and MuToElec().

Referenced by ElecToTau(), Oscillate(), OscillateMatter(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), and Extrapolate2D::SetNominalOscProb().

00459 {
00460   // Flip delta to reverse direction
00461   double oldSinDelta = fSinDCP;
00462   double oldDelta =  fOscPar[OscPar::kDelta];
00463 
00464   fOscPar[OscPar::kDelta] = -oldDelta;
00465   fSinDCP = -oldSinDelta;
00466                                                                                                                        
00467   double prob = OscCalc::MuToElec(x);
00468                                                                                                                        
00469   //restore the world
00470   fOscPar[OscPar::kDelta] = oldDelta;
00471   fSinDCP = oldSinDelta;
00472  
00473   return prob;   
00474 }

double OscCalc::ElecToTau ( double * x )

Definition at line 501 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References ElecToMu(), fCosTh, fSin2Th, and OscPar::kTh23.

00502 {
00503 //  EtoTau is the same as E->Mu wich sinsq_th23 <-> cossq_th23, sin(2th23) <->-sin(2th23)
00504   double origCos = fCosTh[OscPar::kTh23];
00505   double origSin = fSinTh[OscPar::kTh23];
00506   double orig2Sin = fSin2Th[OscPar::kTh23];
00507 
00508   fCosTh[OscPar::kTh23] = origSin;
00509   fSinTh[OscPar::kTh23] = origCos;
00510   fSin2Th[OscPar::kTh23] = -orig2Sin;
00511 
00512   double prob = OscCalc::ElecToMu(x);
00513 
00514   //restore the world
00515   fCosTh[OscPar::kTh23] = origCos;
00516   fSinTh[OscPar::kTh23] = origSin;
00517   fSin2Th[OscPar::kTh23] = orig2Sin;
00518 
00519   return prob;
00520 
00521 /*  This is an option as well but I think results in more cycles
00522   double p1 = 1. - OscCalc::ElecToMu(x) - OscCalc::ElecToElec(x);
00523   if(p1 < 0){ cout<<"Damnation! "<<x[0]<<" "<<p1<<endl; p1 = 0;}
00524   return p1;
00525 */ 
00526 }

double OscCalc::ElecToTau ( double x )

Definition at line 431 of file DataUtil/OscCalc.cxx.

References ElecToMu(), fCosTh, fSin2Th, fSinTh, and OscPar::kTh23.

Referenced by Oscillate(), OscillateMatter(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::SetNominalOscProb(), and TauToElec().

00432 {
00433 //  EtoTau is the same as E->Mu wich sinsq_th23 <-> cossq_th23, sin(2th23) <->-sin(2th23)
00434   double origCos = fCosTh[OscPar::kTh23];
00435   double origSin = fSinTh[OscPar::kTh23];
00436   double orig2Sin = fSin2Th[OscPar::kTh23];
00437 
00438   fCosTh[OscPar::kTh23] = origSin;
00439   fSinTh[OscPar::kTh23] = origCos;
00440   fSin2Th[OscPar::kTh23] = -orig2Sin;
00441 
00442   double prob = OscCalc::ElecToMu(x);
00443 
00444   //restore the world
00445   fCosTh[OscPar::kTh23] = origCos;
00446   fSinTh[OscPar::kTh23] = origSin;
00447   fSin2Th[OscPar::kTh23] = orig2Sin;
00448 
00449   return prob;
00450 
00451 /*  This is an option as well but I think results in more cycles
00452   double p1 = 1. - OscCalc::ElecToMu(x) - OscCalc::ElecToElec(x);
00453   if(p1 < 0){ cout<<"Damnation! "<<x[0]<<" "<<p1<<endl; p1 = 0;}
00454   return p1;
00455 */ 
00456 }

void OscCalc::GetOscParam ( double * par )

double OscCalc::GetOscParam ( OscPar::OscPar_t pos )

Definition at line 98 of file DataUtil/OscCalc.cxx.

References fOscPar, and OscPar::kUnknown.

00099 {
00100    if(pos >=0 && pos < OscPar::kUnknown)
00101      return fOscPar[pos];
00102    else 
00103      return -9999;
00104 }

void OscCalc::GetOscParam ( double * par )

Definition at line 91 of file DataUtil/OscCalc.cxx.

References fOscPar, and OscPar::kUnknown.

Referenced by NueStandard::GetOscParam(), NueSystematic::GetOscParams(), Extrapolate2D::GetPrediction(), Extrapolate2D::InvertMassHierarchy(), Extrapolate2D_Simple::OscillatePrediction(), Extrapolate2D::OscillatePrediction(), and Extrapolate2D::ReadExtrap().

00092 {
00093    for(int i = 0; i < int(OscPar::kUnknown); i++){
00094      par[i] = fOscPar[i];
00095    }
00096 }

double OscCalc::MuSurvive ( double * x )

Definition at line 494 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References MuToElec(), and MuToTau().

Referenced by OscillateMatter().

00495 {
00496   double p1 = 1. - OscCalc::MuToTau(x) - OscCalc::MuToElec(x);
00497   if(p1 < 0){ cout<<"Damnation! "<<x[0]<<" "<<p1<<endl; p1 = 0;}
00498   return p1;
00499 }

double OscCalc::MuToElec ( double * x )

Definition at line 183 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References MuELoss::e, fCos2Th, fCosDCP, fCosTh, fOscPar, fSin2Th, fSinDCP, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, OscPar::kTh13, and OscPar::kTh23.

00184 {
00185   double E = x[0]; //energy
00186 
00187   //Building the standard terms
00188   double L = fOscPar[OscPar::kL];  //baseline
00189   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00190   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00191   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00192 
00193   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00194   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00195                                                                                 
00196   double cos_th23 = fCosTh[OscPar::kTh23];
00197   double cos_th12 = fCosTh[OscPar::kTh12];
00198   double sin_th13 = fSinTh[OscPar::kTh13];
00199                                                                                 
00200   double sin_2th23 = fSin2Th[OscPar::kTh23];
00201   double sin_2th12 = fSin2Th[OscPar::kTh12];
00202   double cos_2th13 = fCos2Th[OscPar::kTh13];
00203   double cos_2th12 = fCos2Th[OscPar::kTh12];
00204 
00205   double sinsq_th23 = fSinTh[OscPar::kTh23]*fSinTh[OscPar::kTh23];
00206   double sinsq_th12 = fSinTh[OscPar::kTh12]*fSinTh[OscPar::kTh12];
00207 
00208 //      printf("josh sin osc par th23 %f val %f\n",fSinTh[OscPar::kTh23],fOscPar[OscPar::kTh23]);
00209 
00210 //      printf("josh sinsq_th23 %f sinsq_2th13 %f\n",sinsq_th23,sinsq_2th13);
00211 
00212   double d_cp = fOscPar[OscPar::kDelta];
00213   double cos_dcp = fCosDCP;
00214   double sin_dcp = fSinDCP;
00215 
00216   //Building the more complicated terms
00217   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00218   double A = 2*fV*E*1e9/(dmsq_13);
00219   double alpha = dmsq_12/dmsq_13;
00220 
00221   // A and d_cp both change sign for antineutrinos
00222   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00223   A *= plusminus;
00224   d_cp *= plusminus;
00225   sin_dcp *= plusminus;
00226 
00227   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00228 
00229   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00230 
00231 //printf("josh C13 = %f has sinsq_2th13 %f cos_2th13 %f A %f\n", C13,sinsq_2th13,cos_2th13,A);
00232 
00233   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00234   if(fabs(alpha) > 1e-10){  //want to be careful here
00235     double temp = cos_2th12 - A/alpha;
00236     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00237   }
00238 
00239   //More complicated sin and cosine terms
00240   double cosC13Delta = cos(C13*Delta);
00241   double sinC13Delta = sin(C13*Delta);
00242   double sin1pADelta = sin((A+1)*Delta);
00243   double cos1pADelta = cos((A+1)*Delta);
00244   double sinADelta = sin(A*Delta);
00245   double sinAm1Delta = sin((A-1)*Delta);
00246   double cosdpDelta = cos(d_cp+Delta);
00247   double sinApam2Delta = sin((A+alpha-2)*Delta);
00248   double cosApam2Delta = cos((A+alpha-2)*Delta);
00249 
00250   double cosaC12Delta = 0;
00251   double sinaC12Delta = 0; 
00252   
00253   if(fabs(alpha) > 1e-10){
00254     cosaC12Delta = cos(alpha*C12*Delta);
00255     sinaC12Delta = sin(alpha*C12*Delta);
00256   } // otherwise not relevant
00257 
00258   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00259   // this is the equivalent of Eq 47 & 48 corrected for Mu to E instead of E to Mu
00260 
00261   // Leading order term 
00262   double p1 = sinsq_th23*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00263 //printf("josh l/e term %f\n",sinC13Delta*sinC13Delta);
00264 //      printf("p1 %f sinsq_th23 %f sinsq_2th13 %f * sinC13Delta^2 %f / C13^2 %f\n",p1,sinsq_th23,sinsq_2th13,sinC13Delta*sinC13Delta,C13*C13);
00265 
00266   // terms that appear at order alpha
00267   double p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 - 
00268                       A*sinC13Delta*(cos_2th13-A)/(C13*C13); 
00269 
00270   double p2 = -2*sinsq_th12*sinsq_th23*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;
00271 
00272   double p3Inner = - sin_dcp*(cosC13Delta - cos1pADelta)*C13 
00273          + cos_dcp*(C13*sin1pADelta - (1-A*cos_2th13)*sinC13Delta);
00274 
00275   double p3 = sin_2th12*sin_2th23*sin_th13*sinC13Delta/(A*C13*C13)*p3Inner*alpha;
00276 
00277   //  p1 + p2 + p3 is the complete contribution for this expansion
00278   
00279   // Now for the expansion in orders of sin(th13) (good to all order alpha) 
00280   //  this is the equivalent of Eq 65 and 66
00281 
00282   // leading order term
00283   double pa1 = 0.0, pa2 = 0.0;
00284 
00285   if(fabs(alpha) > 1e-10){
00286     // leading order term
00287     pa1 = cos_th23*cos_th23*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);
00288 
00289     // and now to calculate the first order in s13 term
00290     double t1 = (cos_2th12 - A/alpha)/C12 
00291                   - alpha*A*C12*sinsq_2th12/(2*(1-alpha)*C12*C12);
00292     double t2 = -cos_dcp*(sinApam2Delta-sinaC12Delta*t1);
00293   
00294     double t3 = -(cosaC12Delta-cosApam2Delta)*sin_dcp;
00295  
00296     double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;
00297     double t4 = sin_2th12*sin_2th23*(1-alpha)*sinaC12Delta/denom;
00298 
00299     pa2 = t4*(t3+t2)*sin_th13;
00300   }
00301   //pa1+pa2 is the complete contribution from this expansion
00302 
00303   // In order to combine the information correctly we need to add the two
00304   //  expansions and subtract off the terms that are in both (alpha^1, s13^1) 
00305   //  these may be taken from the expansion to second order in both parameters
00306   //  Equation 31 
00307 
00308   double t1 = sinADelta*cosdpDelta*sinAm1Delta/(A*(A-1));
00309   double repeated = 2*alpha*sin_2th12*sin_2th23*sin_th13*t1;
00310 
00311   //  Calculate the total probability
00312   double totalP = p1+p2+p3 + (pa1+pa2) - repeated;
00313 
00314 
00315 //      printf("delta %f p1 %f p2 %f p3 %f pa1+pa2 %f repeated %f  prob %f\n",d_cp,p1,p2,p3,pa1+pa2,repeated,totalP);
00316 
00317   return totalP;
00318 }

double OscCalc::MuToElec ( double x )

Definition at line 107 of file DataUtil/OscCalc.cxx.

References MuELoss::e, fCos2Th, fCosDCP, fCosTh, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, OscPar::kTh13, and OscPar::kTh23.

Referenced by ElecToMu(), MuSurvive(), MuToMu(), Oscillate(), OscillateMatter(), Extrapolate2D_Simple::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D_Simple::OscillatePrediction_SepNuNuBar(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), and Extrapolate2D::SetNominalOscProb().

00108 {
00109   double E = x; //energy
00110 
00111   //Building the standard terms
00112   double L = fOscPar[OscPar::kL];  //baseline
00113   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00114   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00115   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00116 
00117   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00118   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00119                                                                                 
00120   double cos_th23 = fCosTh[OscPar::kTh23];
00121   double cos_th12 = fCosTh[OscPar::kTh12];
00122   double sin_th13 = fSinTh[OscPar::kTh13];
00123   double cos_th13 = fCosTh[OscPar::kTh13];  
00124                                                                               
00125   double sin_2th23 = fSin2Th[OscPar::kTh23];
00126   double sin_2th12 = fSin2Th[OscPar::kTh12];
00127   double cos_2th13 = fCos2Th[OscPar::kTh13];
00128   double cos_2th12 = fCos2Th[OscPar::kTh12];
00129 
00130   double sinsq_th23 = fSinTh[OscPar::kTh23]*fSinTh[OscPar::kTh23];
00131   double sinsq_th12 = fSinTh[OscPar::kTh12]*fSinTh[OscPar::kTh12];
00132 
00133   double d_cp = fOscPar[OscPar::kDelta];
00134   double cos_dcp = fCosDCP;
00135   double sin_dcp = fSinDCP;
00136 
00137   //Building the more complicated terms
00138   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00139   double A = 2*fV*E*1e9/(dmsq_13);
00140   double alpha = dmsq_12/dmsq_13;
00141 
00142   // A and d_cp both change sign for antineutrinos
00143   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00144   A *= plusminus;
00145   d_cp *= plusminus;
00146   sin_dcp *= plusminus;
00147 
00148   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00149 
00150   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00151 
00152   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00153   if(fabs(alpha) > 1e-10){  //want to be careful here
00154     double temp = cos_2th12 - A/alpha;
00155     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00156   }
00157 
00158   //More complicated sin and cosine terms
00159   double cosC13Delta = cos(C13*Delta);
00160   double sinC13Delta = sin(C13*Delta);
00161   double sin1pADelta = sin((A+1)*Delta);
00162   double cos1pADelta = cos((A+1)*Delta);
00163   double sinADelta = sin(A*Delta);
00164   double sinAm1Delta = sin((A-1)*Delta);
00165   double cosdpDelta = cos(d_cp+Delta);
00166   double sinApam2Delta = sin((A+alpha-2)*Delta);
00167   double cosApam2Delta = cos((A+alpha-2)*Delta);
00168 
00169   double cosaC12Delta = 0;
00170   double sinaC12Delta = 0; 
00171   
00172   if(fabs(alpha) > 1e-10){
00173     cosaC12Delta = cos(alpha*C12*Delta);
00174     sinaC12Delta = sin(alpha*C12*Delta);
00175   } // otherwise not relevant
00176 
00177   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00178   // this is the equivalent of Eq 47 & 48 corrected for Mu to E instead of E to Mu
00179 
00180   // Leading order term 
00181   double p1 = sinsq_th23*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00182 
00183   // terms that appear at order alpha
00184   //first work out the vacuum case since we get 0/0 otherwise.......
00185   double p2Inner = Delta*cosC13Delta;
00186 
00187   if(fabs(A) > 1e-9)
00188     p2Inner = Delta*cosC13Delta*(1-A*cos_2th13)/C13 -
00189                       A*sinC13Delta*(cos_2th13-A)/(C13*C13);
00190 
00191   double p2 = -2*sinsq_th12*sinsq_th23*sinsq_2th13*sinC13Delta/(C13*C13)*p2Inner*alpha;
00192 
00193 
00194   //again working out vacuum first.....
00195   double p3Inner = Delta* cos_th13* cos_th13*(-2*sin_dcp*sinC13Delta*sinC13Delta+2*cos_dcp*sinC13Delta*cosC13Delta);
00196 
00197   if(fabs(A) > 1e-9)
00198     p3Inner = (sinC13Delta/(A*C13*C13))*(- sin_dcp*(cosC13Delta - cos1pADelta)*C13
00199                          + cos_dcp*(C13*sin1pADelta - (1-A*cos_2th13)*sinC13Delta));
00200 
00201   double p3 = sin_2th12*sin_2th23*sin_th13*p3Inner*alpha;
00202 
00203   //  p1 + p2 + p3 is the complete contribution for this expansion
00204   
00205   // Now for the expansion in orders of sin(th13) (good to all order alpha) 
00206   //  this is the equivalent of Eq 65 and 66
00207 
00208   // leading order term
00209   double pa1 = 0.0, pa2 = 0.0;
00210 
00211   // no problems here when A -> 0
00212   if(fabs(alpha) > 1e-10){
00213     // leading order term
00214     pa1 = cos_th23*cos_th23*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12);
00215 
00216     // and now to calculate the first order in s13 term
00217     double t1 = (cos_2th12 - A/alpha)/C12 
00218                   - alpha*A*C12*sinsq_2th12/(2*(1-alpha)*C12*C12);
00219     double t2 = -cos_dcp*(sinApam2Delta-sinaC12Delta*t1);
00220   
00221     double t3 = -(cosaC12Delta-cosApam2Delta)*sin_dcp;
00222  
00223     double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;
00224     double t4 = sin_2th12*sin_2th23*(1-alpha)*sinaC12Delta/denom;
00225 
00226     pa2 = t4*(t3+t2)*sin_th13;
00227   }
00228   //pa1+pa2 is the complete contribution from this expansion
00229 
00230   // In order to combine the information correctly we need to add the two
00231   //  expansions and subtract off the terms that are in both (alpha^1, s13^1) 
00232   //  these may be taken from the expansion to second order in both parameters
00233   //  Equation 31 
00234 
00235   double t1 = Delta*sinC13Delta*cosdpDelta;
00236   if(fabs(A) > 1e-9) t1 = sinADelta*cosdpDelta*sinAm1Delta/(A*(A-1));
00237 
00238   double repeated = 2*alpha*sin_2th12*sin_2th23*sin_th13*t1;
00239 
00240   //  Calculate the total probability
00241   double totalP = p1+p2+p3 + (pa1+pa2) - repeated;
00242 
00243   return totalP;
00244 }

double OscCalc::MuToMu ( double x )

Definition at line 424 of file DataUtil/OscCalc.cxx.

References MuELoss::e, MuToElec(), and MuToTau().

Referenced by Oscillate(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::SetNominalOscProb(), and TauToTau().

00425 {
00426   double p1 = 1. - OscCalc::MuToTau(x) - OscCalc::MuToElec(x);
00427   if(p1 < 1e-4){ cout<<"P(mu->mu) less than zero Damnation! "<<x<<" "<<p1<<endl; p1 = 0;}
00428   return p1;
00429 }

double OscCalc::MuToTau ( double * x )

Definition at line 320 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References MuELoss::e, fCos2Th, fCosDCP, fCosTh, fOscPar, fSin2Th, fSinDCP, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, OscPar::kTh13, and OscPar::kTh23.

00321 {
00322   double E = x[0]; //energy
00323 
00324   double L = fOscPar[OscPar::kL];  //baseline
00325   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00326   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00327   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00328 
00329   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00330   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00331   double sinsq_2th23 = fSin2Th[OscPar::kTh23]*fSin2Th[OscPar::kTh23];
00332                                                                                 
00333   double cos_th13 = fCosTh[OscPar::kTh13];
00334   double cos_th12 = fCosTh[OscPar::kTh12];
00335   double sin_th12 = fSinTh[OscPar::kTh12];
00336   double sin_th13 = fSinTh[OscPar::kTh13];
00337                                                                                 
00338   double sin_2th23 = fSin2Th[OscPar::kTh23];
00339 //  double sin_2th13 = fSin2Th[OscPar::kTh13];
00340   double sin_2th12 = fSin2Th[OscPar::kTh12];
00341                                                                                 
00342   double cos_2th23 = fCos2Th[OscPar::kTh23];
00343   double cos_2th13 = fCos2Th[OscPar::kTh13];
00344   double cos_2th12 = fCos2Th[OscPar::kTh12];
00345 
00346   double d_cp = fOscPar[OscPar::kDelta];
00347   double cos_dcp = fCosDCP;
00348   double sin_dcp = fSinDCP;
00349 
00350   //Building the more complicated terms                                                                              
00351   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00352   double A = 2*fV*E*1e9/(dmsq_13);
00353   double alpha = dmsq_12/dmsq_13;
00354   
00355   // A and d_cp both change sign for antineutrinos
00356   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00357   A *= plusminus;
00358   d_cp *= plusminus;
00359   sin_dcp *= plusminus;
00360 
00361   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00362                                                                                                                        
00363   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00364                                                                                                                        
00365   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00366   if(fabs(alpha) > 1e-10){  //want to be careful here
00367     double temp = cos_2th12 - A/alpha;
00368     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00369   }
00370                                                                                                                        
00371   //More complicated sin and cosine terms
00372   double cosC13Delta = cos(C13*Delta);
00373   double sinC13Delta = sin(C13*Delta);
00374   double sin1pADelta = sin((A+1)*Delta);
00375   double cos1pADelta = cos((A+1)*Delta);
00376   double sinADelta = sin((A)*Delta);
00377 
00378   double sin1pAmCDelta = sin(0.5*(A+1-C13)*Delta);
00379   double sin1pApCDelta = sin(0.5*(A+1+C13)*Delta);
00380 
00381   double cosaC12Delta = 0;
00382   double sinaC12Delta = 0;
00383                                                                                                                        
00384   if(fabs(alpha) > 1e-10){
00385     cosaC12Delta = cos(alpha*C12*Delta);
00386     sinaC12Delta = sin(alpha*C12*Delta);
00387   } // otherwise not relevant
00388 
00389   double sinApam2Delta = sin((A+alpha-2)*Delta);
00390   double cosApam2Delta = cos((A+alpha-2)*Delta);
00391   double sinAm1Delta = sin((A-1)*Delta);
00392   double cosAm1Delta = cos((A-1)*Delta);
00393   double sinDelta = sin(Delta);
00394   double sin2Delta = sin(2*Delta);
00395 
00396   double cosaC12pApam2Delta = 0;                                                                                                                       
00397   if(fabs(alpha) > 1e-10){
00398     cosaC12pApam2Delta = cos((alpha*C12+A+alpha-2)*Delta);
00399   }
00400 
00401   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00402   // this is the equivalent of Eq 49 & 50 corrected for Mu to E instead of E to Mu
00403 
00404   // Leading order term
00405   double pmt_0 = 0.5*sinsq_2th23;
00406   pmt_0 *= (1 - (cos_2th13-A)/C13)*sin1pAmCDelta*sin1pAmCDelta 
00407              +  (1 + (cos_2th13-A)/C13)*sin1pApCDelta*sin1pApCDelta
00408              - 0.5*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00409 
00410   // terms that appear at order alpha
00411   double t0, t1, t2, t3;
00412   t0 = (cos_th12*cos_th12-sin_th12*sin_th12*sin_th13*sin_th13
00413           *(1+2*sin_th13*sin_th13*A+A*A)/(C13*C13))*cosC13Delta*sin1pADelta*2;
00414   t1 = 2*(cos_th12*cos_th12*cos_th13*cos_th13-cos_th12*cos_th12*sin_th13*sin_th13
00415          +sin_th12*sin_th12*sin_th13*sin_th13
00416          +(sin_th12*sin_th12*sin_th13*sin_th13-cos_th12*cos_th12)*A);
00417   t1 *= sinC13Delta*cos1pADelta/C13;
00418 
00419   t2 =  sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13*C13);
00420   t2 *= A/Delta*sin1pADelta+A/Delta*(cos_2th13-A)/C13*sinC13Delta
00421           - (1-A*cos_2th13)*cosC13Delta;
00422 
00423   double pmt_1 = -0.5*sinsq_2th23*Delta*(t0+t1+t2);   
00424 
00425   t0 = t1 = t2 = t3 = 0.0;
00426 
00427   t0 = cosC13Delta-cos1pADelta;
00428   t1 = 2*cos_th13*cos_th13*sin(d_cp)*sinC13Delta/C13*t0;
00429   t2 = -cos_2th23*cos_dcp*(1+A)*t0*t0;
00430 
00431   t3  = cos_2th23*cos_dcp*(sin1pADelta+(cos_2th13-A)/C13*sinC13Delta);
00432   t3 *= (1+2*sin_th13*sin_th13*A + A*A)*sinC13Delta/C13 - (1+A)*sin1pADelta;
00433 
00434 //  cout<<t1<<"  "<<t2<<"  "<<t3<<endl;
00435 
00436   pmt_1 = pmt_1 + (t1+t2+t3)*sin_th13*sin_2th12*sin_2th23/(2*A*cos_th13*cos_th13);
00437   pmt_1 *= alpha;
00438 
00439   //  pmt_0 + pmt_1 is the complete contribution for this expansion
00440                                                                                                                        
00441   // Now for the expansion in orders of sin(th13) (good to all order alpha)
00442   //  this is the equivalent of Eq 67 and 68
00443                                                                                                                        
00444   // leading order term
00445   double pmt_a0 =  0.5*sinsq_2th23;
00446 
00447   pmt_a0 *= 1 - 0.5*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12)
00448               - cosaC12pApam2Delta
00449               - (1 - (cos_2th12 - A/alpha)/C12)*sinaC12Delta*sinApam2Delta;
00450             
00451   double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;
00452 
00453   t0 = (cosaC12Delta-cosApam2Delta)*(cosaC12Delta-cosApam2Delta);
00454            
00455   t1 = (cos_2th12 - A/alpha)/C12*sinaC12Delta+sinApam2Delta;
00456             
00457   t2 = ((cos_2th12 - A/alpha)/C12+2*(1-alpha)/(alpha*A*C12))*sinaC12Delta
00458              + sinApam2Delta;
00459 
00460   t3 = (alpha*A*C12)/2.0*cos_2th23*cos_dcp*(t0 + t1*t2);
00461   
00462   t3 += sin_dcp*(1-alpha)*(cosaC12Delta-cosApam2Delta)*sinaC12Delta;
00463 
00464   double pmt_a1 = sin_th13*sin_2th12*sin_2th23/denom*t3;
00465 
00466   // pmt_a1+pmt_a2 is the complete contribution from this expansion
00467                                                                                                                        
00468   // In order to combine the information correctly we need to add the two
00469   //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
00470   //  and lower order terms
00471   //  these may be taken from the expansion to second order in both parameters
00472   //  Equation 34
00473 
00474 
00475   // Now for the term of order alpha * s13 or lower order!
00476   t0 = t1 = t2 = t3 = 0.0;
00477 
00478   t1 = +sin_dcp*sinDelta*sinADelta*sinAm1Delta/(A*(A-1));
00479   t2 = -1/(A-1)*cos_dcp*sinDelta*(A*sinDelta-sinADelta*cosAm1Delta/A)*cos_2th23/denom;
00480   t0 =  2*alpha*sin_2th12*sin_2th23*sin_th13*(t1+t2);
00481 
00482   t1 = sinsq_2th23*sinDelta*sinDelta 
00483        - alpha*sinsq_2th23*cos_th12*cos_th12*Delta*sin2Delta;
00484 
00485   double repeated = t0+t1;
00486 
00487   //  Calculate the total probability
00488   double totalP = pmt_0 + pmt_1 + pmt_a0 + pmt_a1 - repeated;
00489 
00490   return totalP;
00491 }

double OscCalc::MuToTau ( double x )

Definition at line 246 of file DataUtil/OscCalc.cxx.

References MuELoss::e, fCos2Th, fCosDCP, fCosTh, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::invKmToeV, OscPar::kDelta, OscPar::kDeltaM12, OscPar::kDeltaM23, OscPar::kL, OscPar::kNuAntiNu, OscPar::kTh12, OscPar::kTh13, and OscPar::kTh23.

Referenced by MuSurvive(), MuToMu(), Oscillate(), OscillateMatter(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::SetNominalOscProb(), and TauToMu().

00247 {
00248   double E = x; //energy
00249 
00250   double L = fOscPar[OscPar::kL];  //baseline
00251   double dmsq_23 = fOscPar[OscPar::kDeltaM23];
00252   double dmsq_12 = fOscPar[OscPar::kDeltaM12];
00253   double dmsq_13 = dmsq_23+dmsq_12; //eV^{2}
00254 
00255   double sinsq_2th12 = fSin2Th[OscPar::kTh12]*fSin2Th[OscPar::kTh12];
00256   double sinsq_2th13 = fSin2Th[OscPar::kTh13]*fSin2Th[OscPar::kTh13];
00257   double sinsq_2th23 = fSin2Th[OscPar::kTh23]*fSin2Th[OscPar::kTh23];
00258                                                                                 
00259   double cos_th13 = fCosTh[OscPar::kTh13];
00260   double cos_th12 = fCosTh[OscPar::kTh12];
00261   double sin_th12 = fSinTh[OscPar::kTh12];
00262   double sin_th13 = fSinTh[OscPar::kTh13];
00263                                                                                 
00264   double sin_2th23 = fSin2Th[OscPar::kTh23];
00265 //  double sin_2th13 = fSin2Th[OscPar::kTh13];
00266   double sin_2th12 = fSin2Th[OscPar::kTh12];
00267                                                                                 
00268   double cos_2th23 = fCos2Th[OscPar::kTh23];
00269   double cos_2th13 = fCos2Th[OscPar::kTh13];
00270   double cos_2th12 = fCos2Th[OscPar::kTh12];
00271 
00272   double d_cp = fOscPar[OscPar::kDelta];
00273   double cos_dcp = fCosDCP;
00274   double sin_dcp = fSinDCP;
00275 
00276   //Building the more complicated terms                                                                              
00277   double Delta = dmsq_13*L/(4*E*1e9*OscPar::invKmToeV);
00278   double A = 2*fV*E*1e9/(dmsq_13);
00279   double alpha = dmsq_12/dmsq_13;
00280   
00281   // A and d_cp both change sign for antineutrinos
00282   double plusminus = int(fOscPar[OscPar::kNuAntiNu]);
00283   A *= plusminus;
00284   d_cp *= plusminus;
00285   sin_dcp *= plusminus;
00286 
00287   //Now calculate the resonance terms for alpha expansion (C13) and s13 expansion (C12)
00288                                                                                                                        
00289   double C13 = TMath::Sqrt(sinsq_2th13+(A-cos_2th13)*(A-cos_2th13));
00290                                                                                                                        
00291   double C12 = 1;  //really C12 -> infinity when alpha = 0 but not an option really
00292   if(fabs(alpha) > 1e-10){  //want to be careful here
00293     double temp = cos_2th12 - A/alpha;
00294     C12 = TMath::Sqrt(sinsq_2th12+(temp*temp));
00295   }
00296                                                                                                                        
00297   //More complicated sin and cosine terms
00298   double cosC13Delta = cos(C13*Delta);
00299   double sinC13Delta = sin(C13*Delta);
00300   double sin1pADelta = sin((A+1)*Delta);
00301   double cos1pADelta = cos((A+1)*Delta);
00302   double sinADelta = sin((A)*Delta);
00303 
00304   double sin1pAmCDelta = sin(0.5*(A+1-C13)*Delta);
00305   double sin1pApCDelta = sin(0.5*(A+1+C13)*Delta);
00306 
00307   double cosaC12Delta = 0;
00308   double sinaC12Delta = 0;
00309                                                                                                                        
00310   if(fabs(alpha) > 1e-10){
00311     cosaC12Delta = cos(alpha*C12*Delta);
00312     sinaC12Delta = sin(alpha*C12*Delta);
00313   } // otherwise not relevant
00314 
00315   double sinApam2Delta = sin((A+alpha-2)*Delta);
00316   double cosApam2Delta = cos((A+alpha-2)*Delta);
00317   double sinAm1Delta = sin((A-1)*Delta);
00318   double cosAm1Delta = cos((A-1)*Delta);
00319   double sinDelta = sin(Delta);
00320   double sin2Delta = sin(2*Delta);
00321 
00322   double cosaC12pApam2Delta = 0;                                                                                                                       
00323   if(fabs(alpha) > 1e-10){
00324     cosaC12pApam2Delta = cos((alpha*C12+A+alpha-2)*Delta);
00325   }
00326 
00327   //First we calculate the terms for the alpha expansion (good to all orders in th13)
00328   // this is the equivalent of Eq 49 & 50 corrected for Mu to E instead of E to Mu
00329 
00330   // Leading order term
00331   double pmt_0 = 0.5*sinsq_2th23;
00332   pmt_0 *= (1 - (cos_2th13-A)/C13)*sin1pAmCDelta*sin1pAmCDelta 
00333              +  (1 + (cos_2th13-A)/C13)*sin1pApCDelta*sin1pApCDelta
00334              - 0.5*sinsq_2th13*sinC13Delta*sinC13Delta/(C13*C13);
00335 
00336   // terms that appear at order alpha
00337   double t0, t1, t2, t3;
00338   t0 = (cos_th12*cos_th12-sin_th12*sin_th12*sin_th13*sin_th13
00339           *(1+2*sin_th13*sin_th13*A+A*A)/(C13*C13))*cosC13Delta*sin1pADelta*2;
00340   t1 = 2*(cos_th12*cos_th12*cos_th13*cos_th13-cos_th12*cos_th12*sin_th13*sin_th13
00341          +sin_th12*sin_th12*sin_th13*sin_th13
00342          +(sin_th12*sin_th12*sin_th13*sin_th13-cos_th12*cos_th12)*A);
00343   t1 *= sinC13Delta*cos1pADelta/C13;
00344 
00345   t2 =  sin_th12*sin_th12*sinsq_2th13*sinC13Delta/(C13*C13*C13);
00346   t2 *= A/Delta*sin1pADelta+A/Delta*(cos_2th13-A)/C13*sinC13Delta
00347           - (1-A*cos_2th13)*cosC13Delta;
00348 
00349   double pmt_1 = -0.5*sinsq_2th23*Delta*(t0+t1+t2);   
00350 
00351   t0 = t1 = t2 = t3 = 0.0;
00352 
00353   t0 = cosC13Delta-cos1pADelta;
00354   t1 = 2*cos_th13*cos_th13*sin(d_cp)*sinC13Delta/C13*t0;
00355   t2 = -cos_2th23*cos_dcp*(1+A)*t0*t0;
00356 
00357   t3  = cos_2th23*cos_dcp*(sin1pADelta+(cos_2th13-A)/C13*sinC13Delta);
00358   t3 *= (1+2*sin_th13*sin_th13*A + A*A)*sinC13Delta/C13 - (1+A)*sin1pADelta;
00359 
00360 //  cout<<t1<<"  "<<t2<<"  "<<t3<<endl;
00361 
00362   if(fabs(A) > 1e-9) 
00363     pmt_1 = pmt_1 + (t1+t2+t3)*sin_th13*sin_2th12*sin_2th23/(2*A*cos_th13*cos_th13);
00364   else
00365     pmt_1 = pmt_1 + sin_th13*sin_2th12*sin_2th23*cos_th13*cos_th13*Delta*(2*sin_dcp*sinC13Delta*sinC13Delta+cos_dcp*cos_2th23*2*sinC13Delta*cosC13Delta);
00366 
00367   pmt_1 *= alpha;
00368 
00369   //  pmt_0 + pmt_1 is the complete contribution for this expansion
00370                                                                                                                        
00371   // Now for the expansion in orders of sin(th13) (good to all order alpha)
00372   //  this is the equivalent of Eq 67 and 68
00373                                                                                                                        
00374   // leading order term
00375   double pmt_a0 =  0.5*sinsq_2th23;
00376 
00377   pmt_a0 *= 1 - 0.5*sinsq_2th12*sinaC12Delta*sinaC12Delta/(C12*C12)
00378               - cosaC12pApam2Delta
00379               - (1 - (cos_2th12 - A/alpha)/C12)*sinaC12Delta*sinApam2Delta;
00380             
00381   double denom = (1-A-alpha+A*alpha*cos_th12*cos_th12)*C12;
00382 
00383   t0 = (cosaC12Delta-cosApam2Delta)*(cosaC12Delta-cosApam2Delta);
00384            
00385   t1 = (cos_2th12 - A/alpha)/C12*sinaC12Delta+sinApam2Delta;
00386             
00387   t2 = ((cos_2th12 - A/alpha)/C12+2*(1-alpha)/(alpha*A*C12))*sinaC12Delta
00388              + sinApam2Delta;
00389 
00390   t3 = (alpha*A*C12)/2.0*cos_2th23*cos_dcp*(t0 + t1*t2);
00391   
00392   t3 += sin_dcp*(1-alpha)*(cosaC12Delta-cosApam2Delta)*sinaC12Delta;
00393 
00394   double pmt_a1 = sin_th13*sin_2th12*sin_2th23/denom*t3;
00395 
00396   // pmt_a1+pmt_a2 is the complete contribution from this expansion
00397                                                                                                                        
00398   // In order to combine the information correctly we need to add the two
00399   //  expansions and subtract off the terms that are in both (alpha^1, s13^1)
00400   //  and lower order terms
00401   //  these may be taken from the expansion to second order in both parameters
00402   //  Equation 34
00403 
00404 
00405   // Now for the term of order alpha * s13 or lower order!
00406   t0 = t1 = t2 = t3 = 0.0;
00407 
00408   t1 = +sin_dcp*sinDelta*sinADelta*sinAm1Delta/(A*(A-1));
00409   t2 = -1/(A-1)*cos_dcp*sinDelta*(A*sinDelta-sinADelta*cosAm1Delta/A)*cos_2th23;
00410   t0 =  2*alpha*sin_2th12*sin_2th23*sin_th13*(t1+t2);
00411 
00412   t1 = sinsq_2th23*sinDelta*sinDelta 
00413        - alpha*sinsq_2th23*cos_th12*cos_th12*Delta*sin2Delta;
00414 
00415   double repeated = t0+t1;
00416 
00417   //  Calculate the total probability
00418   double totalP = pmt_0 + pmt_1 + pmt_a0 + pmt_a1 - repeated;
00419 
00420   return totalP;
00421 }

double OscCalc::Oscillate	(	int	nuFlavor,
		int	nonOscNuFlavor,
		double	Energy
	)

double OscCalc::Oscillate	(	int	nuFlavor,
		int	nonOscNuFlavor,
		float	Energy,
		float	L,
		float	dm2,
		float	theta23,
		float	UE32
	)

Definition at line 62 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References fCosTh, fOscPar, fSin2Th, OscPar::kDeltaM23, OscPar::kL, OscPar::kTh23, and SetOscParam().

00064 {
00065    double par[9] = {0};
00066    for(int i = 0; i < 9; i++) {par[i] = fOscPar[i];}
00067    par[OscPar::kL] = L;
00068    par[OscPar::kDeltaM23] = dm2;
00069    par[OscPar::kTh23] = theta23;
00070 
00071    SetOscParam(par);
00072 
00073    double oscterm = TMath::Sin(1.269*dm2*L/Energy);
00074    oscterm = oscterm*oscterm;
00075 
00076    double sin2th23 = fSin2Th[OscPar::kTh23];
00077    double sinsq2th23 = sin2th23*sin2th23;
00078    double sinth23 = fSinTh[OscPar::kTh23];
00079    double costh23 = fCosTh[OscPar::kTh23];
00080 
00081    double pmt=(1-UE32)*(1-UE32)*sinsq2th23*oscterm;
00082    double pme=sinth23*sinth23*4.*UE32*(1-UE32)*oscterm;
00083    double pmm=1.-pmt-pme;
00084 
00085    double pet=4*(1-UE32)*UE32*oscterm*costh23*costh23;
00086    double pem=sinth23*sinth23*4.*UE32*(1-UE32)*oscterm;
00087    double pee=1.-pet-pem;
00088 
00089 
00090    if(abs(nonOscNuFlavor)==14){
00091       if(abs(nuFlavor)==12)       { return pme; }
00092       else if(abs(nuFlavor)==14)  { return pmm; }
00093       else if(abs(nuFlavor)==16)  { return pmt; }
00094    }
00095    else if(abs(nonOscNuFlavor)==12){
00096       if(abs(nuFlavor)==12)       { return pee; }
00097       else if(abs(nuFlavor)==14)  { return pem; }
00098       else if(abs(nuFlavor)==16)  { return pet; }
00099    }
00100    else{
00101      std::cout<<"I don't know what to do with "<<nonOscNuFlavor
00102          <<" "<<nuFlavor<<" "<<pee<<std::endl;
00103    }
00104    return 0.;
00105 }

double OscCalc::Oscillate	(	int	nuFlavor,
		int	nonOscNuFlavor,
		double	Energy,
		double *	par
	)

Definition at line 19 of file DataUtil/OscCalc.cxx.

References Oscillate(), and SetOscParam().

00020 {
00021   SetOscParam(par);
00022   return Oscillate(nuFlavor, nonOscNuFlavor, Energy); 
00023 }

double OscCalc::Oscillate	(	int	nuFlavor,
		int	nonOscNuFlavor,
		double	Energy
	)

Definition at line 26 of file DataUtil/OscCalc.cxx.

References ElecToElec(), ElecToMu(), ElecToTau(), Msg::kInfo, OscPar::kNuAntiNu, MSG, MuToElec(), MuToMu(), MuToTau(), SetOscParam(), TauToElec(), TauToMu(), and TauToTau().

Referenced by NueSystematic::DoOscCalc(), ANtpTruthInfoBeamAna::GetOscProb(), NueStandard::GetOscWeight(), Oscillate(), NueSensitivity::OscillateMatter(), NueFCSensitivity::OscillateMatter(), NuOscProb::ThreeFlavourNuMuToNuElectron(), NuOscProb::ThreeFlavourNuMuToNuMu(), NuOscProb::ThreeFlavourNuMuToNuTau(), and NC::OscProb::ThreeFlavor::TransitionProbability().

00027 {
00028   double x = Energy;
00029 
00030   if(nonOscNuFlavor > 0) SetOscParam(OscPar::kNuAntiNu, 1);
00031   else SetOscParam(OscPar::kNuAntiNu, -1);
00032 
00033   double prob = 0.0;
00034   if(abs(nonOscNuFlavor)==14) {
00035     if(abs(nuFlavor)==12)      prob = OscCalc::MuToElec(x);
00036     else if(abs(nuFlavor)==14) prob = OscCalc::MuToMu(x);
00037     else if(abs(nuFlavor)==16) prob = OscCalc::MuToTau(x);
00038   }
00039   if(abs(nonOscNuFlavor)==12) {
00040     if(abs(nuFlavor)==14) prob = OscCalc::ElecToMu(x);
00041     else if(abs(nuFlavor)==12) prob = OscCalc::ElecToElec(x);
00042     else if(abs(nuFlavor)==16) prob = OscCalc::ElecToTau(x);
00043   }
00044   if(abs(nonOscNuFlavor)==16) {
00045     if(abs(nuFlavor)==14) prob = OscCalc::TauToMu(x);
00046     else if(abs(nuFlavor)==12) prob = OscCalc::TauToElec(x);
00047     else if(abs(nuFlavor)==16) prob = OscCalc::TauToTau(x);
00048   } 
00049 
00050   if(prob < -1e-4 && Energy > 0.05){  //I'm sorry but the oscillations are just noisy at low E 
00051       MSG("OscCalc",Msg::kInfo)<<"P(<"<<nonOscNuFlavor<<"->"<<nuFlavor
00052                                <<") less than zero at E = "<<x<<" prob = "<<prob<<endl;
00053       prob = 0;
00054   } 
00055 
00056   return prob;
00057 }

double OscCalc::OscillateMatter	(	int	nuFlavor,
		int	nonOscNuFlavor,
		double	Energy
	)

Definition at line 117 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References ElecToElec(), ElecToMu(), ElecToTau(), OscPar::kNuAntiNu, MuSurvive(), MuToElec(), MuToTau(), and SetOscParam().

00119 {
00120   double x[1] = {};
00121   x[0] = Energy;
00122 
00123   if(nonOscNuFlavor > 0) SetOscParam(OscPar::kNuAntiNu, 1);
00124   else SetOscParam(OscPar::kNuAntiNu, -1);
00125 
00126   if(abs(nonOscNuFlavor)==14) {
00127     if(abs(nuFlavor)==12)      return OscCalc::MuToElec(x);
00128     else if(abs(nuFlavor)==14) return OscCalc::MuSurvive(x);
00129     else if(abs(nuFlavor)==16) return OscCalc::MuToTau(x);
00130   }
00131   if(abs(nonOscNuFlavor)==12) {
00132     if(abs(nuFlavor)==14) return OscCalc::ElecToMu(x);
00133     else if(abs(nuFlavor)==12) return OscCalc::ElecToElec(x);
00134     else if(abs(nuFlavor)==16) return OscCalc::ElecToTau(x);
00135   }
00136 
00137   return 0;
00138 }

double OscCalc::OscillateMatter	(	int	nuFlavor,
		int	nonOscNuFlavor,
		double	Energy,
		double *	par
	)

Definition at line 109 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References SetOscParam().

00111 {
00112   SetOscParam(par);
00113   return OscillateMatter(nuFlavor, nonOscNuFlavor, Energy); 
00114 }

void OscCalc::Print ( )

Definition at line 57 of file NueAna/ParticlePID/ParticleAna/OscCalc.cxx.

References fOscPar.

00057                     {
00058         for(int i=0;i<10;i++)printf("%f\n",fOscPar[i]);
00059    
00060    }

void OscCalc::SetOscParam	(	OscPar::OscPar_t	pos,
		double	val,
		bool	force = `false`
	)

void OscCalc::SetOscParam	(	double *	par,
		bool	force = `false`
	)

void OscCalc::SetOscParam	(	OscPar::OscPar_t	pos,
		double	val,
		bool	force = `false`
	)

Definition at line 59 of file DataUtil/OscCalc.cxx.

References OscPar::A_av, MuELoss::e, fCos2Th, fCosDCP, fCosTh, fElecDensity, fOscPar, fSin2Th, fSinDCP, fSinTh, fV, OscPar::Gf, OscPar::invCmToeV, OscPar::invCmToGeV, OscPar::kDelta, OscPar::kDensity, OscPar::root2, and OscPar::z_a.

00060 {
00061   if(fabs(fOscPar[pos] - val) > 1e-9 || force){
00062     fOscPar[pos] = val;
00063     if(pos < 3){
00064        fSinTh[pos] = sin(val);
00065        fCosTh[pos] = cos(val);
00066        fSin2Th[pos] = sin(2*val);
00067        fCos2Th[pos] = cos(2*val);
00068     }
00069     if(pos == OscPar::kDensity){
00070       double ne = OscPar::z_a*OscPar::A_av*val; //electron density #/cm^{3}
00071       fElecDensity = ne*OscPar::invCmToeV*OscPar::invCmToGeV*OscPar::invCmToGeV;
00072       //electron density with units Gev^{2} eV
00073       //Gev^{2} to cancel with GeV^{-2} in Gf
00074                                                                                 
00075        fV = OscPar::root2*OscPar::Gf*fElecDensity; //eV
00076     }
00077     if(pos == OscPar::kDelta){
00078        fCosDCP = cos(val); 
00079        fSinDCP = sin(val);
00080     }
00081   }
00082 }

void OscCalc::SetOscParam	(	double *	par,
		bool	force = `false`
	)

Definition at line 84 of file DataUtil/OscCalc.cxx.

References OscPar::kUnknown.

Referenced by NueSystematic::DoOscCalc(), ANtpTruthInfoBeamAna::GetOscProb(), Extrapolate2D::InitializeOscCalc(), Extrapolate2D::InvertMassHierarchy(), OscCalc(), Oscillate(), OscillateMatter(), NueSensitivity::OscillateMatter(), NueFCSensitivity::OscillateMatter(), Extrapolate2D_Simple::OscillatePrediction_NoSep(), Extrapolate2D::OscillatePrediction_NoSep(), Extrapolate2D_Simple::OscillatePrediction_SepNuNuBar(), Extrapolate2D::OscillatePrediction_SepNuNuBar(), Extrapolate2D::ReadExtrap(), Extrapolate2D_Simple::ReadFNFile(), Extrapolate2D::ReadFNFile(), Extrapolate2D_Simple::ReadFNFileNo3D(), NueSensitivity::RunStandardApproach(), NueFCSensitivity::RunStandardApproach(), NueStandard::SetDefaultOscParam(), Extrapolate2D::SetDeltaCP(), Extrapolate2D::SetNominalOscProb(), NueStandard::SetOscNoMatter(), Extrapolate2D::SetOscPar(), NueStandard::SetOscParam(), NueSensitivity::SetOscParamBase(), NueFCSensitivity::SetOscParamBase(), NNTrain::SetOscParamBase(), NueSystematic::SetOscParams(), Extrapolate2D::SetSinSq2Th13(), NuOscProb::ThreeFlavourNuMuToNuElectron(), NuOscProb::ThreeFlavourNuMuToNuMu(), NuOscProb::ThreeFlavourNuMuToNuTau(), and NC::OscProb::ThreeFlavor::TransitionProbability().

00085 {
00086    for(int i = 0; i < int(OscPar::kUnknown); i++){
00087      SetOscParam(OscPar::OscPar_t(i), par[i], force);
00088    }
00089 }

double OscCalc::TauToElec ( double x )

Definition at line 607 of file DataUtil/OscCalc.cxx.

References ElecToTau(), fOscPar, fSinDCP, and OscPar::kDelta.

Referenced by Oscillate().

00608 {
00609   // Flip delta to reverse direction
00610   double oldSinDelta = fSinDCP;
00611   double oldDelta =  fOscPar[OscPar::kDelta];
00612    
00613   fOscPar[OscPar::kDelta] = -oldDelta;
00614   fSinDCP = -oldSinDelta; 
00615                                
00616   double prob = OscCalc::ElecToTau(x);
00617                                  
00618   //restore the world 
00619   fOscPar[OscPar::kDelta] = oldDelta;
00620   fSinDCP = oldSinDelta; 
00621    
00622   return prob; 
00623 }

double OscCalc::TauToMu ( double x )

Definition at line 589 of file DataUtil/OscCalc.cxx.

References fOscPar, fSinDCP, OscPar::kDelta, and MuToTau().

Referenced by Oscillate().

00590 {
00591   // Flip delta to reverse direction
00592   double oldSinDelta = fSinDCP;
00593   double oldDelta =  fOscPar[OscPar::kDelta];
00594   
00595   fOscPar[OscPar::kDelta] = -oldDelta;
00596   fSinDCP = -oldSinDelta;
00597                               
00598   double prob = OscCalc::MuToTau(x);
00599                                                                                       
00600   //restore the world
00601   fOscPar[OscPar::kDelta] = oldDelta;
00602   fSinDCP = oldSinDelta;
00603   
00604   return prob;
00605 }

double OscCalc::TauToTau ( double x )

Definition at line 568 of file DataUtil/OscCalc.cxx.

References fCosTh, fSin2Th, fSinTh, OscPar::kTh23, and MuToMu().

Referenced by Oscillate().

00569 {
00570 //  TautoTau is the same as Mu->Mu wich sinsq_th23 <-> cossq_th23, sin(2th23) <->-sin(2th23)
00571   double origCos = fCosTh[OscPar::kTh23];
00572   double origSin = fSinTh[OscPar::kTh23];
00573   double orig2Sin = fSin2Th[OscPar::kTh23];
00574 
00575   fCosTh[OscPar::kTh23] = origSin;
00576   fSinTh[OscPar::kTh23] = origCos;
00577   fSin2Th[OscPar::kTh23] = -orig2Sin;
00578 
00579   double prob = OscCalc::MuToMu(x);
00580   
00581   //restore the world
00582   fCosTh[OscPar::kTh23] = origCos;
00583   fSinTh[OscPar::kTh23] = origSin;
00584   fSin2Th[OscPar::kTh23] = orig2Sin;
00585   
00586   return prob;
00587 }

double OscCalc::fCos2Th [private]

Definition at line 81 of file DataUtil/OscCalc.h.

Referenced by ElecToElec(), MuToElec(), MuToTau(), and SetOscParam().

double OscCalc::fCosDCP [private]

Definition at line 85 of file DataUtil/OscCalc.h.

Referenced by MuToElec(), MuToTau(), and SetOscParam().

double OscCalc::fCosTh [private]

Definition at line 79 of file DataUtil/OscCalc.h.

Referenced by ElecToTau(), MuToElec(), MuToTau(), Oscillate(), SetOscParam(), and TauToTau().

double OscCalc::fElecDensity [private]

Definition at line 82 of file DataUtil/OscCalc.h.

Referenced by SetOscParam().

double OscCalc::fOscPar [private]

Definition at line 76 of file DataUtil/OscCalc.h.

Referenced by ElecToElec(), ElecToMu(), GetOscParam(), MuToElec(), MuToTau(), OscCalc(), Oscillate(), Print(), SetOscParam(), TauToElec(), and TauToMu().

double OscCalc::fSin2Th [private]

Definition at line 80 of file DataUtil/OscCalc.h.

Referenced by ElecToElec(), ElecToTau(), MuToElec(), MuToTau(), Oscillate(), SetOscParam(), and TauToTau().

double OscCalc::fSinAL [private]

Definition at line 44 of file NueAna/ParticlePID/ParticleAna/OscCalc.h.

double OscCalc::fSinDCP [private]

Definition at line 84 of file DataUtil/OscCalc.h.

Referenced by ElecToElec(), ElecToMu(), MuToElec(), MuToTau(), SetOscParam(), TauToElec(), and TauToMu().

double OscCalc::fSinTh [private]

Definition at line 78 of file DataUtil/OscCalc.h.

Referenced by ElecToElec(), ElecToTau(), MuToElec(), MuToTau(), SetOscParam(), and TauToTau().

double OscCalc::fV [private]

Definition at line 83 of file DataUtil/OscCalc.h.

Referenced by ElecToElec(), MuToElec(), MuToTau(), and SetOscParam().

OscCalc Class Reference

Public Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Member Functions
	OscCalc ()
void	GetOscParam (double *par)
double	GetOscParam (OscPar::OscPar_t pos)
void	SetOscParam (double *par, bool force=false)
void	SetOscParam (OscPar::OscPar_t pos, double val, bool force=false)
double	Oscillate (int nuFlavor, int nonOscNuFlavor, double Energy)
double	Oscillate (int nuFlavor, int nonOscNuFlavor, double Energy, double *par)
double	MuToElec (double x)
double	MuToTau (double x)
double	MuToMu (double x)
double	ElecToElec (double x)
double	ElecToMu (double x)
double	ElecToTau (double x)
double	TauToElec (double x)
double	TauToTau (double x)
double	TauToMu (double x)
	OscCalc ()
void	GetOscParam (double *par)
void	SetOscParam (double *par, bool force=false)
void	SetOscParam (OscPar::OscPar_t pos, double val, bool force=false)
double	Oscillate (int nuFlavor, int nonOscNuFlavor, float Energy, float L, float dm2, float theta23, float UE32)
double	Oscillate (int nuFlavor, int nonOscNuFlavor, double Energy)
double	OscillateMatter (int nuFlavor, int nonOscNuFlavor, double Energy, double *par)
double	OscillateMatter (int nuFlavor, int nonOscNuFlavor, double Energy)
double	MuToElec (double *x)
double	MuToTau (double *x)
double	MuSurvive (double *x)
double	ElecToElec (double *x)
double	ElecToMu (double *x)
double	ElecToTau (double *x)
void	Print ()
Private Attributes
double	fOscPar [10]
double	fSinTh [3]
double	fCosTh [3]
double	fSin2Th [3]
double	fCos2Th [3]
double	fElecDensity
double	fV
double	fSinDCP
double	fCosDCP
double	fSinAL