CalcEM.cc File Reference

#include "TF2.h"
#include "TROOT.h"
#include "TMath.h"
#include "rte.cc"
#include "KeyFunc.cc"
#include <cmath>
#include <iostream>
#include <map>
#include <vector>
#include "BinFluctuationEM.h"

Go to the source code of this file.

Defines
#define	D_EFF 6.0
#define	X0_EFF 4.1542
#define	RM_EFF 4.0717
#define	ST_WID 4.1
Functions
std::vector< rte >	RTCalc (double, double, double, double, double, double &)
std::map< int, spef >	PSCalc (std::vector< rte >, int, double, double, double, double, double, double, double)
Double_t *	Rotation3D (double , double )
Double_t	Shw2DProfSam (double , double )
Double_t	LongShwProfSam (double , double )
Double_t	TransShwProfSam (double , double )
map< int, spef >	PSCalc (vector< rte > rtmap, int vtxview, double stripUvtx, double stripVvtx, double planevtx, double angle_u, double angle_v, double delta_Phi, double sumE)

Define Documentation

#define D_EFF 6.0

Definition at line 14 of file CalcEM.cc.

Referenced by AlgSubShowerSR::CalculateEnergyVertexAngle(), LongShwProfSamAng(), PSCalc(), StpFluctuaionAng(), AlgSubShowerSR::SubShowerID(), and TranShwProfSamAng().

#define RM_EFF 4.0717

Definition at line 16 of file CalcEM.cc.

Referenced by AlgSubShowerSR::CalculateEnergyVertexAngle(), LongShwProfSamAng(), PSCalc(), and TranShwProfSamAng().

#define ST_WID 4.1

Definition at line 17 of file CalcEM.cc.

Referenced by AlgSubShowerSR::CalculateEnergyVertexAngle(), and PSCalc().

#define X0_EFF 4.1542

Definition at line 15 of file CalcEM.cc.

Referenced by AlgSubShowerSR::CalculateEnergyVertexAngle(), LongShwProfSamAng(), PSCalc(), and TranShwProfSamAng().

Function Documentation

Double_t LongShwProfSam	(	double *	x,
		double *	par
	)

Definition at line 85 of file CalcEM.cc.

References MinosMaterial::Z().

Referenced by Shw2DProfSam().

00085                                               {
00086 
00087   Double_t t = x[0];
00088   Double_t E = par[0]*1000.; //Energy in MeV
00089   if(E<=0) return 0;
00090 
00091   Double_t Z = 24.9814;
00092   Double_t Ec = 20.7941;
00093   Double_t Fs = 0.692366;
00094   Double_t ehat = 0.8752;
00095   
00096   Double_t lny = TMath::Log(E/Ec);
00097   Double_t alpha = 0.21 + (0.492+2.38/Z)*lny;
00098   Double_t T = lny - 0.858;
00099 
00100   T = T - 0.59/Fs - 0.53*(1-ehat);
00101   alpha = alpha - 0.444/Fs;
00102 
00103   Double_t belta  = (alpha-1)/T;
00104   Double_t f = TMath::Power((belta*t),(alpha-1))
00105     *belta*TMath::Exp(-(belta*t))/TMath::Gamma(alpha);
00106 
00107   return f;
00108 
00109 }

map<int,spef> PSCalc	(	vector< rte >	rtmap,
		int	vtxview,
		double	stripUvtx,
		double	stripVvtx,
		double	planevtx,
		double	angle_u,
		double	angle_v,
		double	delta_Phi,
		double	sumE
	)

Definition at line 160 of file CalcEM.cc.

References D_EFF, MakeKey(), RM_EFF, Rotation3D(), ST_WID, and X0_EFF.

00164 00165 00166 00167 00168 00169 00170 00171 00172 00173 00174 00175 00176 00177 00178 00179 00180 00181 00182 00183 00184 00185 00186 00187 00188 00189 00190 00191 00192 00193 00194 00195 00196 00197 00198 00199 00200 00201 00202 00203 00204 00205       } 00206 00207 00208 00209       } 00210 00211 00212     } 00213     beg++; 00214   } 00215 00216 }

class="fragment">00163 { vector<rte>::iterator beg = rtmap.begin(); vector<rte>::iterator end = rtmap.end(); double pars[2]; pars[0] = TMath::ATan(TMath::Tan(angle_u)*(ST_WID*X0_EFF)/(D_EFF*RM_EFF)); pars[1] = TMath::ATan(TMath::Tan(angle_v)*(ST_WID*X0_EFF)/(D_EFF*RM_EFF)); map<int,spef> psvals; while(beg!=end){ double rtpcoord[3] = {0}; rtpcoord[0] = beg->t; rtpcoord[1] = beg->r; for(double phi=0;phi<2.*TMath::Pi();phi+=delta_Phi){ rtpcoord[2] = phi; double *uvzcoord = Rotation3D(rtpcoord,pars); int signcoord0 = 1; double stU_d = stripUvtx + uvzcoord[0]*RM_EFF/ST_WID; if(stU_d<0) signcoord0 = -1; int stU = int(stU_d + signcoord0*0.5); int signcoord1 = 1; double stV_d = stripVvtx + uvzcoord[1]*RM_EFF/ST_WID; if(stV_d<0) signcoord1 = -1; int stV = int(stV_d + signcoord1*0.5); int signcoord2 = 1; double pl_d = planevtx + uvzcoord[2]*X0_EFF/D_EFF; if(pl_d<0) signcoord2 = -1; int pl = int(pl_d + signcoord2*0.5); delete [] uvzcoord; int key = 0; if(pl%2==0) { //this plane has same plane-view as vtx if(vtxview==2) key = MakeKey(pl,stU); //U-strips else if(vtxview==3) key = MakeKey(pl,stV); //V-strips else { //this plane has diff plane-view to vtx if(vtxview==2) key = MakeKey(pl,stV); else if(vtxview==3) key = MakeKey(pl,stU); psvals[key].en += beg->e*delta_Phi/(2.*TMath::Pi()*sumE); psvals[key].fl += beg->f*delta_Phi/(2.*TMath::Pi()*sumE); return psvals;

std::map<int,spef> PSCalc	(	std::vector< rte >	,
		int	,
		double	,
		double	,
		double	,
		double	,
		double	,
		double	,
		double
	)

Referenced by FitterEM::PredictEMLoss().

Double_t * Rotation3D	(	double *	x,
		double *	par
	)

Definition at line 219 of file CalcEM.cc.

Referenced by PSCalc().

00219                                             {
00220 
00221   Double_t t = x[0];
00222   Double_t r = x[1];
00223   Double_t phi = x[2];
00224   Double_t theta_u = par[0];
00225   Double_t theta_v = par[1];
00226   
00227   Double_t u1,v1,z1;
00228   Double_t u2,v2,z2;
00229   Double_t u,v,z;
00230 
00231   u1 = r*cos(phi);
00232   v1 = r*sin(phi);
00233   z1 = t;
00234 
00235   u2 = u1;
00236   v2 = v1*cos(-theta_v)-z1*sin(-theta_v);
00237   z2 = v1*sin(-theta_v)+z1*cos(-theta_v);
00238 
00239   u = u2*cos(-theta_u)-z2*sin(-theta_u);
00240   v = v2;
00241   z = u2*sin(-theta_u)+z2*cos(-theta_u);
00242 
00243   Double_t *uvz = new Double_t[3];
00244   uvz[0] = u;
00245   uvz[1] = v;
00246   uvz[2] = z;
00247   return uvz;
00248 }

vector< rte > RTCalc	(	double	energy,
		double	delta_T,
		double	delta_R,
		double	epsilon,
		double	cutoff,
		double &	sumE
	)

Definition at line 30 of file CalcEM.cc.

References BinFluctuationEM::CalcFluctuation(), rte::e, rte::f, rte::r, MuELoss::R, Shw2DProfSam(), and rte::t.

Referenced by FitterEM::PredictEMLoss().

00031                                                              {
00032   
00033   BinFluctuationEM *binEM = new BinFluctuationEM(energy);
00034   
00035   TF2 *prof = NULL;
00036   if(gROOT->FindObject("FITSHOWEREM_2DPROFILE")) {
00037     prof = (TF2*) gROOT->FindObject("FITSHOWEREM_2DPROFILE");
00038   }
00039   else {
00040     prof = new TF2("FITSHOWEREM_2DPROFILE",Shw2DProfSam,0.,20.,0.,20.,1);
00041   }
00042   prof->SetParameter(0,energy);
00043   
00044   vector<rte> rtvals;
00045   bool keepGoingInT = true;
00046 
00047   sumE = 0;
00048   double T = 0;
00049   double R = 0;
00050   int counter = 0;
00051   while(keepGoingInT){
00052     bool keepGoingInR = true;
00053     R=0;
00054     while(keepGoingInR){
00055       double fracE = prof->Integral(T,T+delta_T,R,R+delta_R,epsilon);
00056       if(fracE<cutoff||fracE!=fracE) {
00057         keepGoingInR = false;
00058         continue;
00059       }
00060       rte vals;
00061       vals.t = T+delta_T/2.;
00062       vals.r = R+delta_R/2.;
00063       vals.e = fracE;
00064       vals.f = vals.e*binEM->CalcFluctuation(vals.t,vals.r);
00065       sumE+=fracE;
00066       rtvals.push_back(vals);
00067       R+=delta_R;
00068       counter++;
00069     }
00070     if(R==0) keepGoingInT = false;
00071     T+=delta_T;
00072   }
00073   delete binEM;
00074   return rtvals;
00075 }

Double_t Shw2DProfSam	(	double *	x,
		double *	par
	)

Definition at line 78 of file CalcEM.cc.

References LongShwProfSam(), and TransShwProfSam().

Referenced by RTCalc().

00078                                             {
00079   
00080   return LongShwProfSam(x,par)*TransShwProfSam(x,par);
00081   
00082 }

Double_t TransShwProfSam	(	double *	x,
		double *	par
	)

Definition at line 112 of file CalcEM.cc.

References MinosMaterial::Z().

Referenced by Shw2DProfSam().

00112                                                 {
00113 
00114   Double_t t = x[0];
00115   Double_t r = fabs(x[1]);
00116   Double_t E = par[0]*1000.; //Energy in MeV
00117   if(E<=0) return 0;
00118 
00119   Double_t Z = 24.9814;
00120   Double_t Ec = 21.6297;
00121   Double_t Fs = 0.692366;
00122   Double_t ehat = 0.8752;
00123   
00124   Double_t lny = log(E/Ec);
00125   Double_t T = lny - 0.858;
00126   T = T - 0.59/Fs - 0.53*(1-ehat);
00127   Double_t tau = t/T;
00128 
00129   Double_t z1 = 0.0251+0.00319*log(E);
00130   Double_t z2 = 0.1162+(-0.000381*Z);
00131 
00132   Double_t k1 = 0.659+(-0.00309*Z);
00133   Double_t k2 = 0.645;
00134   Double_t k3 = -2.59;
00135   Double_t k4 = 0.3585+0.0421*log(E);
00136 
00137   Double_t p1 = 2.632+(-0.00094*Z);
00138   Double_t p2 = 0.401+0.00187*Z;
00139   Double_t p3 = 1.313+(-0.0686*log(E));
00140 
00141   Double_t Rc = z1+z2*tau;
00142   Double_t Rt = k1*(exp(k3*(tau-k2))+exp(k4*(tau-k2)));
00143   Double_t p = p1*exp((p2-tau)/p3-exp((p2-tau)/p3));
00144 
00145   Rc = Rc - 0.0203*(1.-ehat)+(0.0397/Fs)*exp(-tau);
00146   Rt = Rt - 0.14*(1-ehat)-(0.495/Fs)*exp(-tau);
00147   p = p + (1-ehat)*(0.348-(0.642/Fs)*exp(-TMath::Power(tau-1,2)));
00148 
00149   Double_t fc = 0;
00150   Double_t ft = 0;
00151   if(r!=0) {
00152     fc = 2*r*Rc*Rc/TMath::Power(r*r+Rc*Rc,2);
00153     ft = 2*r*Rt*Rt/TMath::Power(r*r+Rt*Rt,2);
00154   }
00155   Double_t f = p*fc+(1-p)*ft;
00156   return f;
00157 }

CalcEM.cc File Reference

Defines

Functions

Define Documentation

Function Documentation