Skip to content

Commit 4df07bf

Browse files
author
Francesca Di Lodovico
committed
first version
0 parents  commit 4df07bf

17 files changed

+3012
-0
lines changed

‎BargerPropagator.cc

+334
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,334 @@
1+
#include "BargerPropagator.h"
2+
3+
BargerPropagator::BargerPropagator()
4+
{
5+
Earth = new EarthDensity( );
6+
init();
7+
}
8+
9+
10+
BargerPropagator::BargerPropagator( bool k )
11+
{
12+
Earth = new EarthDensity( );
13+
init();
14+
}
15+
16+
17+
BargerPropagator::~BargerPropagator( )
18+
{
19+
delete Earth;
20+
}
21+
22+
BargerPropagator::BargerPropagator( const char * f )
23+
{
24+
Earth = new EarthDensity( f );
25+
init();
26+
}
27+
28+
void BargerPropagator::init()
29+
{
30+
kUseMassEigenstates = false;
31+
32+
//rad earth in [cm] /
33+
REarth = Earth->GetEarthRadiuskm() * 1.0e5;
34+
ProductionHeight = 0.0;
35+
PathLength = 0.0;
36+
37+
// default is neutral matter
38+
density_convert = 0.5;
39+
40+
kAntiMNSMatrix = false ;
41+
kSuppressWarnings = false ;
42+
43+
kOneDominantMass = true ;
44+
}
45+
46+
47+
48+
void BargerPropagator::propagate( int NuFlavor ){
49+
50+
int i,j;
51+
int Layers;
52+
double TransitionMatrix[3][3][2];
53+
double TransitionProduct[3][3][2];
54+
double TransitionTemp[3][3][2];
55+
double RawInputPsi[3][2];
56+
double OutputPsi[3][2];
57+
58+
59+
if( ! kSuppressWarnings )
60+
if(
61+
( kAntiMNSMatrix && NuFlavor > 0) ||
62+
(!kAntiMNSMatrix && NuFlavor < 0)
63+
)
64+
{
65+
std::cout << " Warning BargerPropagator::propagate - " << std::endl;
66+
std::cout << " Propagating neutrino flavor and MNS matrix definition differ :" << std::endl;
67+
std::cout << " MNS Matrix was defined for : " << ( kAntiMNSMatrix ? " Nubar " : "Nu" )<< std::endl;
68+
std::cout << " Propagation is for : " << ( NuFlavor < 0 ? " Nubar " : "Nu" )<< std::endl;
69+
std::cout << " Please check your call to BargerPropagator::SetMNS() " << std::endl;
70+
std::cout << " This message can be suppressed with a call to BargerPropagator::SuppressWarnings() " << std::endl;
71+
72+
exit(-1);
73+
}
74+
75+
clear_complex_matrix( TransitionMatrix );
76+
clear_complex_matrix( TransitionProduct );
77+
clear_complex_matrix( TransitionTemp );
78+
79+
ClearProbabilities();
80+
81+
82+
Earth->SetDensityProfile( CosineZenith, PathLength, ProductionHeight );
83+
Layers = Earth->get_LayersTraversed( );
84+
85+
86+
for ( i = 0; i < Layers ; i++ )
87+
{
88+
get_transition_matrix( NuFlavor,
89+
Energy , // in GeV
90+
Earth->get_DensityInLayer(i) * density_convert,
91+
Earth->get_DistanceAcrossLayer(i)/1.0e5, // in km
92+
TransitionMatrix, // Output transition matrix
93+
0.0 // phase offset
94+
);
95+
96+
if ( i == 0 )
97+
copy_complex_matrix( TransitionMatrix , TransitionProduct );
98+
99+
if ( i >0 ){
100+
clear_complex_matrix( TransitionTemp );
101+
multiply_complex_matrix( TransitionMatrix, TransitionProduct, TransitionTemp );
102+
copy_complex_matrix( TransitionTemp, TransitionProduct );
103+
}//for other layers
104+
}// end of layer loop
105+
106+
107+
// loop on neutrino types
108+
for ( i = 0 ; i < 3 ; i++ )
109+
{
110+
for ( j = 0 ; j < 3 ; j++ )
111+
{ RawInputPsi[j][0] = 0.0; RawInputPsi[j][1] = 0.0; }
112+
113+
if( kUseMassEigenstates )
114+
convert_from_mass_eigenstate( i+1, NuFlavor, RawInputPsi );
115+
else
116+
RawInputPsi[i][0] = 1.0;
117+
118+
119+
multiply_complex_matvec( TransitionProduct, RawInputPsi, OutputPsi );
120+
Probability[i][0] += OutputPsi[0][0] * OutputPsi[0][0] + OutputPsi[0][1]*OutputPsi[0][1];
121+
Probability[i][1] += OutputPsi[1][0] * OutputPsi[1][0] + OutputPsi[1][1]*OutputPsi[1][1];
122+
Probability[i][2] += OutputPsi[2][0] * OutputPsi[2][0] + OutputPsi[2][1]*OutputPsi[2][1];
123+
124+
}//end of neutrino loop
125+
126+
127+
}
128+
129+
130+
131+
132+
void BargerPropagator::ClearProbabilities()
133+
{
134+
for ( int i = 0 ; i < 3; i++ )
135+
for ( int j = 0 ; j < 3 ; j++ )
136+
Probability[i][j] = 0.0;
137+
138+
}
139+
140+
void BargerPropagator::SetMNS( double x12, double x13, double x23,
141+
double m21, double mAtm, double delta,
142+
double Energy_ , bool kSquared, int kNuType )
143+
{
144+
Energy = Energy_;
145+
146+
double sin12;
147+
double sin13;
148+
double sin23;
149+
150+
double lm32 = mAtm ;
151+
// Dominant Mixing mode assumes the user
152+
// simply changes the sign of the input atmospheric
153+
// mixing to invert the hierarchy
154+
// so the input for NH corresponds to m32
155+
// and the input for IH corresponds to m31
156+
if( kOneDominantMass )
157+
{
158+
// For the inverted Hierarchy, adjust the input
159+
// by the solar mixing (should be positive)
160+
// to feed the core libraries the correct value of m32
161+
if( mAtm < 0.0 )
162+
lm32 = mAtm - m21 ;
163+
}
164+
else
165+
{
166+
if( !kSuppressWarnings )
167+
{
168+
std::cout << " BargerPropagator::SetMNS - " << std::endl;
169+
std::cout << " You have opted to specify the value of m23 by yourself. " << std::endl;
170+
std::cout << " This means you must correct the value of m23 when switching " << std::endl;
171+
std::cout << " between the mass hierarchy options. " << std::endl;
172+
std::cout << " This message can be suppressed with BargerPropagator::SuppressWarnings()"<< std::endl;
173+
}
174+
}
175+
176+
177+
178+
//if xAB = sin( xAB )^2
179+
if ( kSquared ){
180+
sin12 = sqrt( x12 );
181+
sin13 = sqrt( x13 );
182+
sin23 = sqrt( x23 );
183+
}
184+
else
185+
{
186+
//if xAB = sin( 2 xAB )^2
187+
sin12 = sqrt( 0.5*(1 - sqrt(1 - x12 )) );
188+
sin13 = sqrt( 0.5*(1 - sqrt(1 - x13 )) );
189+
sin23 = sqrt( 0.5*(1 - sqrt(1 - x23 )) );
190+
}
191+
192+
if ( kNuType < 0 )
193+
{
194+
delta *= -1.0 ;
195+
kAntiMNSMatrix = true ;
196+
}
197+
else
198+
{
199+
kAntiMNSMatrix = false ;
200+
}
201+
202+
init_mixing_matrix( m21, lm32, sin12, sin23, sin13, delta );
203+
204+
}
205+
206+
void BargerPropagator::DefinePath(double cz, double ProdHeight, bool kSetProfile )
207+
{
208+
209+
ProductionHeight = ProdHeight*1e5;
210+
PathLength = sqrt( (REarth + ProductionHeight )*(REarth + ProductionHeight)
211+
- (REarth*REarth)*( 1 - cz*cz)) - REarth*cz;
212+
CosineZenith = cz;
213+
if( kSetProfile )
214+
Earth->SetDensityProfile( CosineZenith, PathLength, ProductionHeight );
215+
216+
}
217+
218+
219+
void BargerPropagator::SetMatterPathLength()
220+
{
221+
222+
int Layers = Earth->get_LayersTraversed( );
223+
224+
MatterPathLength = 0.0;
225+
AirPathLength = 0.0;
226+
for( int i = 1 ; i < Layers ; i++ )
227+
MatterPathLength += Earth->get_DistanceAcrossLayer(i);
228+
229+
AirPathLength += Earth->get_DistanceAcrossLayer(0);
230+
}
231+
232+
void BargerPropagator::SetAirPathLength(double x)
233+
{
234+
// argument is [km], convert to [cm]
235+
AirPathLength = x*1.0e5 - MatterPathLength;
236+
}
237+
238+
239+
240+
double BargerPropagator::GetVacuumProb( int Alpha, int Beta , double Energy, double Path )
241+
{
242+
// alpha -> 1:e 2:mu 3:tau
243+
// Energy[GeV]
244+
// Path[km]
245+
/// simple referes to the fact that in the 3 flavor analysis
246+
// the solar mass term is zero
247+
double Probs[3][3];
248+
249+
250+
get_vacuum_probability( Alpha, Energy, Path, Probs );
251+
252+
Alpha = abs(Alpha);
253+
Beta = abs(Beta);
254+
255+
if ( Alpha > 0 )
256+
return Probs[Alpha-1][Beta-1];
257+
258+
if ( Alpha < 0 ) // assuming CPT!!!
259+
return Probs[Beta-1][Alpha-1];
260+
261+
std::cerr << " BargerPropagator::GetVacuumProb neutrino must be non-zero: " << std::endl;
262+
return -1.0;
263+
264+
}
265+
266+
267+
268+
void BargerPropagator::propagateLinear( int NuFlavor, double pathlength, double Density )
269+
{
270+
271+
int i,j;
272+
273+
double TransitionMatrix[3][3][2];
274+
double TransitionProduct[3][3][2];
275+
double TransitionTemp[3][3][2];
276+
double RawInputPsi[3][2];
277+
double OutputPsi[3][2];
278+
279+
if( ! kSuppressWarnings )
280+
if(
281+
( kAntiMNSMatrix && NuFlavor > 0) ||
282+
(!kAntiMNSMatrix && NuFlavor < 0)
283+
)
284+
{
285+
std::cout << " Warning BargerPropagator::propagateLinear - " << std::endl;
286+
std::cout << " Propagating neutrino flavor and MNS matrix definition differ :" << std::endl;
287+
std::cout << " MNS Matrix was defined for : " << ( kAntiMNSMatrix ? " Nubar " : "Nu" )<< std::endl;
288+
std::cout << " Propagation is for : " << ( NuFlavor < 0 ? " Nubar " : "Nu" )<< std::endl;
289+
std::cout << " Please check your call to BargerPropagator::SetMNS() " << std::endl;
290+
std::cout << " This message can be suppressed with a call to BargerPropagator::SuppressWarnings() " << std::endl;
291+
292+
exit(-1);
293+
}
294+
295+
clear_complex_matrix( TransitionMatrix );
296+
clear_complex_matrix( TransitionProduct );
297+
clear_complex_matrix( TransitionTemp );
298+
299+
ClearProbabilities();
300+
301+
302+
get_transition_matrix( NuFlavor,
303+
Energy , // in GeV
304+
Density * density_convert,
305+
pathlength , // in km
306+
TransitionMatrix, // Output transition matrix
307+
0.0
308+
);
309+
310+
copy_complex_matrix( TransitionMatrix , TransitionProduct );
311+
312+
for ( i = 0 ; i < 3 ; i++ )
313+
{
314+
for ( j = 0 ; j < 3 ; j++ )
315+
{ RawInputPsi[j][0] = 0.0; RawInputPsi[j][1] = 0.0; }
316+
317+
if( kUseMassEigenstates )
318+
convert_from_mass_eigenstate( i+1, NuFlavor, RawInputPsi );
319+
else
320+
RawInputPsi[i][0] = 1.0;
321+
322+
multiply_complex_matvec( TransitionProduct, RawInputPsi, OutputPsi );
323+
324+
Probability[i][0] += OutputPsi[0][0] * OutputPsi[0][0] + OutputPsi[0][1]*OutputPsi[0][1];
325+
Probability[i][1] += OutputPsi[1][0] * OutputPsi[1][0] + OutputPsi[1][1]*OutputPsi[1][1];
326+
Probability[i][2] += OutputPsi[2][0] * OutputPsi[2][0] + OutputPsi[2][1]*OutputPsi[2][1];
327+
328+
}// end of loop on neutrino types
329+
330+
}
331+
332+
333+
334+

0 commit comments

Comments
 (0)