loon: CalcEM.cc Source File

00001 #ifndef CALCEM_C
00002 #define CALCEM_C
00003 #include "TF2.h"
00004 #include "TROOT.h"
00005 #include "TMath.h"
00006 #include "rte.cc"
00007 #include "KeyFunc.cc"
00008 #include <cmath>
00009 #include <iostream>
00010 #include <map>
00011 #include <vector>
00012 #include "BinFluctuationEM.h"
00013 
00014 #define D_EFF 6.0
00015 #define X0_EFF 4.1542
00016 #define RM_EFF 4.0717
00017 #define ST_WID 4.1
00018 
00019 std::vector<rte> RTCalc(double,double,double,double,double,double &);
00020 std::map<int,spef> PSCalc(std::vector<rte>,int,double,double,double,
00021                           double,double,double,double);
00022 Double_t *Rotation3D(double *, double *);
00023 Double_t Shw2DProfSam(double *,double *);
00024 Double_t LongShwProfSam(double *,double *);
00025 Double_t TransShwProfSam(double *,double *);
00026 
00027 using namespace std;
00028 
00029 //*************************************************
00030 vector<rte> RTCalc(double energy,double delta_T,double delta_R,
00031                    double epsilon,double cutoff,double &sumE){
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 }
00076 
00077 //*************************************************
00078 Double_t Shw2DProfSam(double *x,double *par){
00079   
00080   return LongShwProfSam(x,par)*TransShwProfSam(x,par);
00081   
00082 }
00083 
00084 //*************************************************
00085 Double_t LongShwProfSam(double *x,double *par){
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 }
00110 
00111 //**************************************************
00112 Double_t TransShwProfSam(double *x, double *par){
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 }
00158 
00159 //*************************************************
00160 map<int,spef> PSCalc(vector<rte> rtmap,int vtxview,
00161                      double stripUvtx,double stripVvtx,double planevtx,
00162                      double angle_u,double angle_v,double delta_Phi,
00163                      double sumE){
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 }
00217 
00218 //*************************************************
00219 Double_t *Rotation3D(double *x, double *par){
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 }
00249 
00250 #endif //CALCEM_C