66
77DetectorK fat;
88void diagonalise (lutEntry_t &lutEntry);
9+ static float etaMaxBarrel = 1.75 ;
910
1011bool
1112fatSolve (float *eff, float *covm, float pt = 0.1 , float eta = 0.0 , float mass = 0.13957000 , int layer = 0 , int what = 0 , int efftype = 0 , int q = 1 )
@@ -48,8 +49,59 @@ fwdSolve(float *covm, float pt = 0.1, float eta = 0.0, float mass = 0.13957000)
4849 return true ;
4950}
5051
52+ bool
53+ fwdPara (float *covm, float pt = 0.1 , float eta = 0.0 , float mass = 0.13957000 , float Bfield = 0.5 , int layer = 0 , int what = 0 )
54+ {
55+ // parametrised forward response; interpolates between FAT at eta = 1.75 and a fixed parametrisation at eta = 4; only diagonal elements
56+ float eff = 0 ;
57+ float covmbarrel[15 ] = {0 };
58+ if (fabs (eta) < etaMaxBarrel || fabs (eta) > 4 )
59+ return false ;
60+
61+ if (!fatSolve (&eff, covmbarrel, pt, etaMaxBarrel, mass, layer, what)) return false ;
62+
63+ // parametrisation at eta = 4
64+ double beta = 1 ./sqrt (1 +mass*mass/pt/pt/cosh (eta)/cosh (eta));
65+ float dca_pos = 2.5e-4 /sqrt (3 ); // 2.5 micron/sqrt(3)
66+ float r0 = 0.5 ; // layer 0 radius [cm]
67+ float r1 = 1.3 ;
68+ float r2 = 2.5 ;
69+ float x0layer = 0.001 ; // material budget (rad length) per layer
70+ double sigma_alpha = 0.0136 /beta/pt*sqrt (x0layer*cosh (eta))*(1 +0.038 *log (x0layer*cosh (eta)));
71+ double dcaxy_ms = sigma_alpha*r0*sqrt (1 +r1*r1/(r2-r0)/(r2-r0));
72+ double dcaxy2 = dca_pos*dca_pos+dcaxy_ms*dcaxy_ms;
73+
74+ double dcaz_ms = sigma_alpha*r0*cosh (eta);
75+ double dcaz2 = dca_pos*dca_pos+dcaz_ms*dcaz_ms;
76+
77+ float Leta = 2.8 /sinh (eta)-0.01 *r0; // m
78+ double relmomres_pos = 10e-6 *pt/0.3 /Bfield/Leta/Leta*sqrt (720 ./15 .);
79+
80+ float relmomres_barrel = sqrt (covmbarrel[14 ])*pt;
81+ float Router = 1 ; // m
82+ float relmomres_pos_barrel = 10e-6 *pt/0.3 /Bfield/Router/Router/sqrt (720 ./15 .);
83+ float relmomres_MS_barrel = sqrt (relmomres_barrel*relmomres_barrel-relmomres_pos_barrel*relmomres_pos_barrel);
84+
85+ // interpolate MS contrib (rel resolution 0.4 at eta = 4)
86+ float relmomres_MS = 1 ./beta*0.4 *pow (0.4 /(beta*relmomres_MS_barrel),(fabs (eta)-4 .)/(4 .-etaMaxBarrel));
87+ float momres_tot = pt*sqrt (relmomres_pos*relmomres_pos + relmomres_MS*relmomres_MS); // total absolute mom reso
88+
89+ // Fill cov matrix diag
90+ for (int i = 0 ; i < 15 ; ++i)
91+ covm[i] = 0 ;
92+
93+ covm[0 ] = covmbarrel[0 ];
94+ if (dcaxy2 > covm[0 ]) covm[0 ] = dcaxy2;
95+ covm[2 ] = covmbarrel[2 ];
96+ if (dcaz2 > covm[2 ]) covm[2 ] = dcaz2;
97+ covm[5 ] = covmbarrel[5 ]; // sigma^2 sin(phi)
98+ covm[9 ] = covmbarrel[9 ]; // sigma^2 tanl
99+ covm[14 ] = momres_tot*momres_tot/pt/pt/pt/pt; // sigma^2 1/pt
100+ return true ;
101+ }
102+
51103void
52- lutWrite (const char *filename = " lutCovm.dat" , int pdg = 211 , float field = 0.2 , int layer = 0 , int what = 0 , int efftype = 0 )
104+ lutWrite (const char *filename = " lutCovm.dat" , int pdg = 211 , float field = 0.2 , int layer = 0 , int what = 0 , int efftype = 0 , int usePara = 0 )
53105{
54106
55107 // output file
@@ -78,9 +130,9 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2,
78130 lutHeader.radmap .max = 100 .;
79131 // eta
80132 lutHeader.etamap .log = false ;
81- lutHeader.etamap .nbins = 40 ;
82- lutHeader.etamap .min = -2 .;
83- lutHeader.etamap .max = 2 .;
133+ lutHeader.etamap .nbins = 80 ;
134+ lutHeader.etamap .min = -4 .;
135+ lutHeader.etamap .max = 4 .;
84136 // pt
85137 lutHeader.ptmap .log = true ;
86138 lutHeader.ptmap .nbins = 200 ;
@@ -108,7 +160,7 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2,
108160 for (int ipt = 0 ; ipt < npt; ++ipt) {
109161 lutEntry.pt = lutHeader.ptmap .eval (ipt);
110162 lutEntry.valid = true ;
111- if (fabs (eta) < 2 . ) {
163+ if (fabs (eta) <= etaMaxBarrel ) { // full lever arm ends at etaMaxBarrel
112164 // printf(" --- fatSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field);
113165 if (!fatSolve (&lutEntry.eff , lutEntry.covm , lutEntry.pt , lutEntry.eta , lutHeader.mass , layer, what, efftype, q)) {
114166 // printf(" --- fatSolve: error \n");
@@ -120,11 +172,18 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2,
120172 else {
121173 // printf(" --- fwdSolve: pt = %f, eta = %f, mass = %f, field=%f \n", lutEntry.pt, lutEntry.eta, lutHeader.mass, lutHeader.field);
122174 lutEntry.eff = 1 .;
123- if (!fwdSolve (lutEntry.covm , lutEntry.pt , lutEntry.eta , lutHeader.mass )) {
124- // printf(" --- fwdSolve: error \n");
125- lutEntry.valid = false ;
126- for (int i = 0 ; i < 15 ; ++i)
127- lutEntry.covm [i] = 0 .;
175+ bool retval = true ;
176+ if (usePara) {
177+ retval = fwdPara (lutEntry.covm , lutEntry.pt , lutEntry.eta , lutHeader.mass , field, layer, what);
178+ }
179+ else {
180+ retval = fwdSolve (lutEntry.covm , lutEntry.pt , lutEntry.eta , lutHeader.mass );
181+ }
182+ if ( !retval) {
183+ // printf(" --- fwdSolve: error \n");
184+ lutEntry.valid = false ;
185+ for (int i = 0 ; i < 15 ; ++i)
186+ lutEntry.covm [i] = 0 .;
128187 }
129188 }
130189 diagonalise (lutEntry);
0 commit comments