ROCsoft: beamtime/cern-oct12/go4/AnalysisMacros/PR_tools.cxx (r4864/r4738) Source File

Go to the documentation of this file.
00001 #include "TMath.h"
00002 
00003 
00004 Double_t calcSimpleDisplacement(Double_t W, Double_t Q[3]){
00005   Double_t d = -1.0E6;
00006   if (Q[0] > 0.0 && Q[2] > 0.0 && Q[1] > 0.0){
00007     d = 0.5 * W * TMath::Log(Q[2] / Q[0]) / TMath::Log((Q[1] * Q[1]) / (Q[0] * Q[2]));
00008   }
00009   return d;
00010 }
00011 
00012 Double_t calcDisplacement(Double_t W, Double_t Q[3], Double_t sigma){
00013   Double_t d = -1.0E6;
00014   if (Q[0] > 0.0 && Q[2] > 0.0 && Q[1] > 0.0
00015       //&& Q[0] != Q[2]
00016       ){
00017     d = 1. / ((Q[0] * Q[0]) + (Q[2] * Q[2])) * (
00018                                                 (Q[0] * Q[0]) *
00019                                                 ((sigma * sigma) / W * TMath::Log(Q[1] / Q[0]) - (0.5 * W)) +
00020                                                 (Q[2] * Q[2]) *
00021                                                 ((sigma * sigma) / W * TMath::Log(Q[2] / Q[1]) + (0.5 * W))
00022                                                 );
00023   }
00024   return d;
00025 }
00026 Double_t calcDisplacementOverflow(Double_t W, Double_t Q[3]){
00027   Double_t d = -1.0E6;
00028   Double_t av = 1775;//V
00029   Double_t a = -0.0012 * av + 0.14;//[-1.0, -1.3] -> [950V, 1200V]  1775V in 2012  => f(V) = -0.0012/V * V + 0.14  => f(1775V) = -1.99
00030   a = -1.1; 
00031   d = (0.5 * W + a * TMath::Log(Q[0] / Q[2])) - 0.5 * W;
00032   //a = -1.4;
00033   //d = (0.5 * W + a * TMath::Log10(Q[0] / Q[2])) - 0.5 * W;
00034   
00035   return d;
00036 }
00037 Double_t calcDisplacementRatio(Double_t W, Double_t Q[3]){
00038   Double_t d = -1.0E6;
00039   if (Q[0] > 0.0 && Q[2] > 0.0 && Q[1] > 0.0
00040       //&& Q[0] != Q[2]
00041       ){
00042     d = 0.5 * W * (Q[2] - Q[0]) / (Q[1] - Q[0]);
00043   }
00044   return d;
00045 }
00046 Double_t calcDisplacementGRatio(Double_t W, Double_t Q[3]){
00047   Double_t f_l = Q[0] / Q[1];
00048   Double_t f_r = Q[2] / Q[1];
00049   Double_t alpha = TMath::ATan((1. - f_l) / (1. - f_r));
00050   Double_t d = 0.5 * W * TMath::ATan(1.65 * (alpha - 0.25 * TMath::Pi())) / TMath::ATan(1.65 * 0.25 * TMath::Pi());  
00051   return d;
00052 }
00053 
00054 Double_t calcDisplacementHyperbolicSecantSquared(Double_t W, Double_t Q[3]){
00055   Double_t a3 = TMath::Pi() * W / (TMath::ACosH(0.5 * (TMath::Sqrt(Q[1] / Q[0]) + TMath::Sqrt(Q[1] / Q[2]))));
00056   Double_t d = a3 / TMath::Pi() * TMath::ATanH((TMath::Sqrt(Q[1] / Q[0]) - TMath::Sqrt(Q[1] / Q[2])) / (2. * TMath::SinH(TMath::Pi() * W / a3))); 
00057   return d;
00058 }
00059 Double_t calcDisplacementHyperbolicSecantSquaredFitting(Double_t W, Double_t Q[3]){
00060   Double_t a3 = TMath::Pi() * W / (TMath::ACosH(0.5 * (TMath::Sqrt(Q[1] / Q[0]) + TMath::Sqrt(Q[1] / Q[2]))));
00061   TH1D signal ("signal","signal",3,-1.5*W,1.5*W);
00062   //signal.Reset();
00063   for (Int_t i = 0; i < 3; i++)
00064     signal.Fill(-W+i*W,Q[i]);
00065   TString formula;
00066   formula.Form("[0] / (TMath::CosH(TMath::Pi() * (x-[1]) / %E) * TMath::CosH(TMath::Pi() * (x-[1]) / %E))",a3,a3);
00067   TF1 f("f",formula,-1.5*W,1.5*W);
00068   f.SetParameter(0,Q[1]);
00069   signal.Fit("f","R0Q");
00070   Double_t d = f.GetParameter(1);
00071  
00072   return d;
00073 }
00074 Double_t calcDisplacementMathiesonFitting(Double_t W, Double_t Q[3]){
00075   Double_t K3(0.35), h(3.5), d(0.0/*calcSimpleDisplacement(W, Q)*/), Qsum(Q[0] + Q[1] + Q[2]);//, c(0.5 * TMath::Pi() * (1. - 0.5 * TMath::Sqrt(K3))), c2(c / h), c3(c * TMath::Sqrt(K3) / (4. * TMath::ATan(TMath::Sqrt(K3))));
00076   //if (fabs(d) > 0.20 * W)
00077     {
00078       TH1D signal ("signal","signal",3,-1.5*W,1.5*W);
00079       //signal.Reset();
00080       for (Int_t i = 0; i < 3; i++)
00081         signal.Fill( -W + i * W, Q[i] / (Qsum));
00082       TString formula;
00083       //formula.Form("[0] * (1. - (TMath::TanH([1] * x / %E) * TMath::TanH([1] * x / %E))) / (1. + [2] * (TMath::TanH([1] * x / %E) * TMath::TanH([1] * x / %E)))",h,h,h,h);
00084       //formula.Form("[0] * (1. - (TMath::TanH([1] * x / %E) * TMath::TanH([1] * x / %E))) / (1. + %E * (TMath::TanH([1] * x / %E) * TMath::TanH([1] * x / %E)))",h,h,K3,h,h);
00085       //formula.Form("[0] * %E * (1. - (TMath::TanH((x + [1]) * %E) * TMath::TanH((x + [1]) * %E))) / (1. + %E * (TMath::TanH((x + [1]) * %E) * TMath::TanH((x + [1]) * %E)))",c3,c2,c2,K3,c2,c2); // (x+[1])???
00086       formula.Form("[0] * (0.5 * TMath::Pi() *(1. - 0.5 * TMath::Sqrt([2])) * TMath::Sqrt([2])) / (4. * TMath::ATan(TMath::Sqrt([2]))) * (1. - TMath::Power(TMath::TanH(0.5 * TMath::Pi() * (1. - 0.5 * TMath::Sqrt([2])) * (x - [1]) / %E),2)) / (1. + [2] * TMath::Power(TMath::TanH(0.5 * TMath::Pi() * (1. - 0.5 * TMath::Sqrt([2])) * (x - [1]) / %E),2))",h,h);
00087       TF1 f("f",formula,-1.5*W,1.5*W);
00088       f.SetParameter(0,Q[1] / (Qsum));
00089       f.SetParameter(1,d);
00090       f.SetParameter(2,K3);
00091       f.SetParLimits(0, 0.85 * Q[1] / (Qsum), 1.15 * Q[1] / (Qsum));
00092       f.SetParLimits(1,-0.51*W,0.51*W);
00093       f.SetParLimits(2,0.1,0.9);
00094       signal.Fit("f","R0QL");
00095       d = f.GetParameter(1);
00096     }
00097   return d;
00098 }
00099 Double_t calcDisplacementPRFFitting(Double_t W, Double_t Q[3]){
00100   Double_t h(3.5), K3(0.35), Qsum(Q[0] + Q[1] + Q[2]);//, c(TMath::ATan(TMath::Sqrt(K3))), sqrtK3(TMath::Sqrt(K3)), c2(TMath::Pi() * (sqrtK3 - 2.) / 8.);
00101   TH1D signal ("signal","signal",3,-1.5*W,1.5*W);
00102   //signal.Reset();
00103   for (Int_t i = 0; i < 3; i++)
00104     signal.Fill(-W+i*W, Q[i] / (Qsum));
00105   TString formula;
00106   //formula.Form("[0] * (-0.5 * (TMath::ATan([1] * TMath::TanH([2] * %E - 2. * [3] * x / %E)))) + (-0.5 * (TMath::ATan([1] * TMath::TanH([2] * %E + 2. * [3] * x / %E))))",W,h,W,h);
00107   //formula.Form("%E * (-0.5 * (TMath::ATan([0] * TMath::TanH([1] * %E - 2. * [2] * x / %E)))) + (-0.5 * (TMath::ATan([0] * TMath::TanH([1] * %E + 2. * [2] * x / %E))))",c,W,h,W,h);
00108   //formula.Form("[1] * %E * (-0.5 * (TMath::ATan(%E * TMath::TanH(%E * %E - 2. * (x + [0]) / %E)))) + (-0.5 * (TMath::ATan(%E * TMath::TanH(%E * %E + 2. * (x + [0]) / %E))))",c,sqrtK3,c2,W,h,sqrtK3,c2,W,h);// (x+[1])???
00109   formula.Form("1. / (TMath::ATan(TMath::Sqrt([1]))) * ((-0.5 * (TMath::ATan(TMath::Sqrt([1]) * TMath::TanH(TMath::Pi() * (TMath::Sqrt([1]) - 2.) * (%E - 2. * (x - [0])) / (8. * %E))))) + (-0.5 * (TMath::ATan(TMath::Sqrt([1]) * TMath::TanH(TMath::Pi() * (TMath::Sqrt([1]) - 2.) * (%E + 2. * (x - [0])) / (8. * %E))))))",W,h,W,h);
00110   TF1 f("f",formula,-1.5*W,1.5*W);
00111   //f.SetParameter(0,Q[1]/Qsum);
00112   f.SetParameter(1,K3);
00113   f.SetParLimits(1,0.1,0.9);
00114 
00115   f.SetParameter(0,0.0);
00116   f.SetParLimits(0,-0.51*W,0.51*W);
00117   signal.Fit("f","R0Q");
00118   //Double_t d = f.GetParameter(3);
00119   //Double_t d = f.GetParameter(2);
00120   Double_t d = f.GetParameter(0);
00121   return d;
00122 }
00123 
00124 void PR_tools(){
00125 
00126 }