StRoot: StEEmcPool/StEEmcPointMaker/eeTowerFunction.cxx Source File

00001 #include "eeTowerFunction.h"
00002 #include "TSpline.h"
00003 #include "TMath.h" 
00004 Double_t eeTowerFunction( Double_t *x, Double_t *p ) {
00005 
00007 
00008   Double_t p0[]={ 0.235772, 0.239713, 0.239243, 0.238843, 0.24141,  0.241641, 0.2433,   0.243124, 0.242957, 0.243987, 0.244477, 0.245425, 0.245425};
00009   Double_t p2[]={ 0.266629, 0.283157, 0.270271, 0.257208, 0.246998, 0.236122, 0.226747, 0.218976, 0.21176,  0.205661, 0.198209, 0.197987, 0.197987};
00010   Double_t p3[]={ 0.497183, 0.490622, 0.491357, 0.491976, 0.492601, 0.493278, 0.493775, 0.494394, 0.495469, 0.494036, 0.493174, 0.504341, 0.504341};
00011 
00012   Double_t xx[]={ 0.,1.,2.,3.,4.,5.,6.,7.,8.,9.,10.,11.,12.};
00013 
00014   /*
00015   TSpline3 *mEtaPar0=new TSpline3("EtaPar0",xx,p0,13);
00016   TSpline5 *mEtaPar2=new TSpline5("EtaPar2",xx,p2,13);
00017   TSpline5 *mEtaPar3=new TSpline5("EtaPar3",xx,p3,13);
00018   */
00019   TSpline3 mEtaPar0("EtaPar0",xx,p0,13);
00020   TSpline5 mEtaPar2("EtaPar2",xx,p2,13);
00021   TSpline5 mEtaPar3("EtaPar3",xx,p3,13);
00022 
00023   Double_t pp0[]={ 0.23439,  0.242071, 0.243334, 0.243623, 0.245357, 0.245606, 0.245668, 0.247,    0.247328, 0.247872, 0.248534, 0.247648, 0.247648 };
00024   Double_t pp2[]={ 0.296042, 0.273496, 0.248509, 0.227491, 0.208265, 0.195849, 0.189432, 0.176265, 0.166762, 0.151045, 0.14633,  0.136524, 0.136524 };
00025   Double_t pp3[]={ 0.490751, 0.490351, 0.492535, 0.495746, 0.496635, 0.498088, 0.497785, 0.498563, 0.499,    0.499639, 0.499539, 0.499887, 0.499887 };
00026   TSpline3 mPhiPar0("PhiPar0",xx,pp0,13);
00027   TSpline3 mPhiPar2("PhiPar2",xx,pp2,13);
00028   TSpline3 mPhiPar3("PhiPar3",xx,pp3,13);
00029 
00030   Double_t xphi=x[0];
00031   Double_t xeta=x[1];
00032 
00033   Double_t xphi0 = p[0];
00034   Double_t xeta0 = p[1];
00035   Double_t energy = p[2];
00036 
00037 
00038   Double_t norm = energy * ( mPhiPar0 .Eval(xeta) * mEtaPar0 .Eval(xeta) );
00039 
00040   Double_t phi_par2 = mPhiPar2 .Eval( xeta );
00041   Double_t eta_par2 = mEtaPar2 .Eval( xeta );
00042 
00043   Double_t phi_par3 = mPhiPar3 .Eval( xeta );
00044   Double_t eta_par3 = 0.50; // Fixed, since it varies so little.
00045 
00046   //  Double_t norm = energy * 0.25;
00047   //  Double_t phi_par2=0.27;
00048   //uble_t eta_par2=0.30;
00049   //  Double_t phi_par3=0.50;
00050   //  Double_t eta_par3=0.49;
00051 
00053   Double_t phi1 = TMath::Erfc( -(xphi-xphi0+phi_par3) / phi_par2 );
00054   Double_t phi2 = TMath::Erfc( +(xphi-xphi0-phi_par3) / phi_par2 );
00055 
00057   Double_t eta1 = TMath::Erfc( -(xeta-xeta0+eta_par3) / eta_par2 );
00058   Double_t eta2 = TMath::Erfc( +(xeta-xeta0-eta_par3) / eta_par2 );
00059 
00061   Double_t result = norm * phi1*phi2 * eta1*eta2;
00062 
00063   return result;
00064 
00065 };
00066 
00067 
00068 
00069 
00070 Double_t eeTower2Function( Double_t *x, Double_t *p ) {
00071 
00072 
00073   Double_t phibin=p[0]; // Fractional phi of higher E gamma
00074   Double_t etabin=p[1]; // Fractional eta of higher E gamma
00075   Double_t energy=p[2]; // Energy of higher E gamma
00076 
00077   Double_t dgg=p[3];    // gamma-gamma separation in eta-phi space
00078   Double_t theta=p[4];  // rotation angle in eta-phi space
00079   Double_t zgg=p[5];    // energy sharing between gammas
00080 
00081   Double_t phibin2=phibin+TMath::Cos(theta) * dgg;
00082   Double_t etabin2=etabin+TMath::Sin(theta) * dgg;
00083   Double_t energy2=energy * (1.0-zgg)/(1.0+zgg);
00084 
00085   Double_t p1[]={phibin,etabin,energy};
00086   Double_t p2[]={phibin2,etabin2,energy2};
00087 
00088   return eeTowerFunction(x,p1) + eeTowerFunction(x,p2);
00089 
00090 };