(*For faster calculation, assume n=5 scintillators (in the experiment there were 10)*)
n=5;
(*R1 =((Scintillator response (T1, 1 photon))/X), where X = M (measured data) or δM (measurement error)*)
R1=Array[r1,n]
(*The same for T2*)
R2=Array[r2,n]
(*The same for T3*)
R3=Array[r3,n]
(*Q=(M/X); NB: Q = {1,1,1,1,1} for X = M*)
Q=Array[q,n]

(*For three temperatures*)

dif=Q-n1 R1-n2 R2-n3 R3
mf3=dif.dif

d1=(1/2)Simplify[(∂_n1) mf3]
d2=(1/2)Simplify[(∂_n2) mf3]
d3=(1/2)Simplify[(∂_n3) mf3]

Collect[d1,{n1,n2,n3}]
Print[]
Collect[d2,{n1,n2,n3}]
Print[]
Collect[d3,{n1,n2,n3}]

(*Check*)
Simplify[d1==-Q.R1+n1 R1.R1 + n2 R1.R2 + n3 R1.R3]
Simplify[d2==-Q.R2+n1 R2.R1 + n2 R2.R2 + n3 R2.R3]
Simplify[d3==-Q.R3+n1 R3.R1 + n2 R3.R2 + n3 R3.R3]

s3=FullSimplify[Solve[{
QR1==n1 R1R1 + n2 R1R2 + n3 R1R3,
QR2==n1 R1R2 + n2 R2R2 + n3 R2R3,
QR3==n1 R1R3 + n2 R2R3 + n3 R3R3},
{n1,n2,n3}]][[1]]

(*Check the solution with the substitututions*)
subst={QR1->Q.R1,QR2->Q.R2,QR3->Q.R3,R1R1->R1.R1,R1R2->R1.R2,R1R3->R1.R3,R2R2->R2.R2,R2R3->R2.R3,R3R3->R3.R3};
Simplify[d1/.s3/.subst]
Simplify[d2/.s3/.subst]
Simplify[d3/.s3/.subst]

(*For c/c++ codes*)
CForm[n1/.s3]

CForm[n2/.s3]

CForm[n3/.s3]

(*Check that the denominators in the equations for n1, n2, n3 are equal*)
Simplify[Denominator[n1/.s3]==Denominator[n2/.s3]]

Simplify[Denominator[n1/.s3]==Denominator[n3/.s3]]

(*For two temperatures*)

dif=Q-n1 R1-n2 R2
mf2=dif.dif

d1=(1/2)Simplify[(∂_n1) mf2]
d2=(1/2)Simplify[(∂_n2) mf2]

Collect[d1,{n1,n2}]
Print[]
Collect[d2,{n1,n2}]

(*Check*)
Simplify[d1==-Q.R1+n1 R1.R1 + n2 R1.R2]
Simplify[d2==-Q.R2+n1 R2.R1 + n2 R2.R2 ]

s2=FullSimplify[Solve[{
QR1==n1 R1R1 + n2 R1R2 ,
QR2==n1 R1R2 + n2 R2R2 },
{n1,n2}]][[1]]

(*Check the solution with the substitututions*)
subst={QR1->Q.R1,QR2->Q.R2,R1R1->R1.R1,R1R2->R1.R2,R2R2->R2.R2};
Simplify[d1/.s2/.subst]
Simplify[d2/.s2/.subst]

CForm[n1/.s2]

CForm[n2/.s2]

(*Check that the denominators in the equations for n1, n2, n3 are equal*)
Denominator[n1/.s2]==-Denominator[n2/.s2]

(*For one temperature*)

dif=Q-n1 R1
mf1=dif.dif

d1=(1/2)Simplify[(∂_n1) mf1]

Collect[d1,{n1,n2}]

Simplify[d1==-Q.R1+n1 R1.R1 ]

s1=FullSimplify[Solve[{
QR1==n1 R1R1 },
{n1}]][[1]]

(*Check the solution with the substitututions*)
subst={QR1->Q.R1,R1R1->R1.R1};
Simplify[d1/.s1/.subst]

CForm[n1/.s1]