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< rteRTCalc (double, double, double, double, double, double &)
std::map< int, spefPSCalc (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, spefPSCalc (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
#define RM_EFF   4.0717
#define ST_WID   4.1

Definition at line 17 of file CalcEM.cc.

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

#define X0_EFF   4.1542

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.

00163                                  {
00164 
00165   vector<rte>::iterator beg = rtmap.begin();
00166   vector<rte>::iterator end = rtmap.end();
00167   
00168   double pars[2];
00169   pars[0] = TMath::ATan(TMath::Tan(angle_u)*(ST_WID*X0_EFF)/(D_EFF*RM_EFF));
00170   pars[1] = TMath::ATan(TMath::Tan(angle_v)*(ST_WID*X0_EFF)/(D_EFF*RM_EFF));
00171 
00172   map<int,spef> psvals;
00173 
00174   while(beg!=end){
00175     
00176     double rtpcoord[3] = {0};
00177     rtpcoord[0] = beg->t;
00178     rtpcoord[1] = beg->r;
00179     
00180     for(double phi=0;phi<2.*TMath::Pi();phi+=delta_Phi){
00181       rtpcoord[2] = phi;
00182       double *uvzcoord = Rotation3D(rtpcoord,pars);
00183 
00184       int signcoord0 = 1;
00185       double stU_d = stripUvtx + uvzcoord[0]*RM_EFF/ST_WID;
00186       if(stU_d<0) signcoord0 = -1;
00187       int stU = int(stU_d + signcoord0*0.5);
00188 
00189       int signcoord1 = 1;
00190       double stV_d = stripVvtx + uvzcoord[1]*RM_EFF/ST_WID;
00191       if(stV_d<0) signcoord1 = -1;
00192       int stV = int(stV_d + signcoord1*0.5);
00193 
00194       int signcoord2 = 1;
00195       double pl_d = planevtx + uvzcoord[2]*X0_EFF/D_EFF;
00196       if(pl_d<0) signcoord2 = -1;
00197       int pl = int(pl_d + signcoord2*0.5);
00198       
00199       delete [] uvzcoord;
00200       
00201       int key = 0;
00202       if(pl%2==0) {  //this plane has same plane-view as vtx
00203         if(vtxview==2)       key = MakeKey(pl,stU); //U-strips      
00204         else if(vtxview==3)  key = MakeKey(pl,stV); //V-strips
00205       }
00206       else {         //this plane has diff plane-view to vtx
00207         if(vtxview==2)       key = MakeKey(pl,stV);
00208         else if(vtxview==3)  key = MakeKey(pl,stU);
00209       }
00210       psvals[key].en += beg->e*delta_Phi/(2.*TMath::Pi()*sumE);
00211       psvals[key].fl += beg->f*delta_Phi/(2.*TMath::Pi()*sumE);
00212     }
00213     beg++;
00214   }
00215   return psvals;
00216 }

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 }


Generated on 24 Apr 2018 for loon by  doxygen 1.6.1