| 1 |
#include <TrkNuclei.h> |
| 2 |
|
| 3 |
|
| 4 |
//////////////////////////////////////////////////////////// |
| 5 |
/// functions |
| 6 |
//////////////////////////////////////////////////////////// |
| 7 |
|
| 8 |
/* |
| 9 |
* First charge-vs-beta parameterization provided by wolfgang |
| 10 |
*/ |
| 11 |
Double_t charge_tracker_beta_wolfgang(Double_t *x, Double_t *par){ |
| 12 |
|
| 13 |
Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 }; |
| 14 |
Float_t B0[9] = {0 , 5.21108 , 8.53345 , 7.00253 , 8.34315 , 11.9997 , 11.6236 , 2.59707 , 2000 }; |
| 15 |
Float_t B1[9] = {0 , -9.6716 , 9.79003 , 38.0199 , 42.6475 , 29.6548 , 43.2637 , 107.077 , 0 }; |
| 16 |
Float_t B2[9] = {0 , 5.70349 , -36.0582 , -73.3067 , -69.9404 , -38.9489 , -60.9272 , -152.77 , 0 }; |
| 17 |
Float_t B3[9] = {0 , -0.311956 , 21.7421 , 37.4468 , 33.7097 , 17.7431 , 31.2389 , 75.0796 , 0 }; |
| 18 |
Float_t x1[9],y1[9]; |
| 19 |
Int_t n1 = 9; |
| 20 |
|
| 21 |
Float_t charge = 1000.; |
| 22 |
Float_t beta = (Float_t) x[0]; |
| 23 |
Float_t dedxtrk = (Float_t) x[1]; dedxtrk = dedxtrk/1.3; |
| 24 |
|
| 25 |
if( beta>2. || dedxtrk <=0 )return 0.;// it makes no sense to allow beta>2 |
| 26 |
|
| 27 |
Float_t betahelp = pow(beta, 0.8); |
| 28 |
Float_t ym = dedxtrk*betahelp; |
| 29 |
Float_t xb = beta; |
| 30 |
for ( Int_t jj=0; jj<9; jj++ ){ |
| 31 |
x1[jj] = B0[jj]+B1[jj]*xb+B2[jj]*xb*xb+B3[jj]*xb*xb*xb; |
| 32 |
y1[jj] = C0[jj]*C0[jj] ; |
| 33 |
} |
| 34 |
TGraph *gr1 = new TGraph(n1,x1,y1); |
| 35 |
TSpline3 *spl1 = new TSpline3("grs",gr1); // use a cubic spline |
| 36 |
Float_t chelp = spl1->Eval(ym); |
| 37 |
charge = TMath::Sqrt(chelp); |
| 38 |
|
| 39 |
if(gr1)gr1->Delete(); |
| 40 |
if(spl1)spl1->Delete(); |
| 41 |
|
| 42 |
return (Double_t)charge ; |
| 43 |
|
| 44 |
} |
| 45 |
/* |
| 46 |
* First charge-vs-rigidity parameterization provided by wolfgang |
| 47 |
* (adpted from charge-vs-deflection parameterization) |
| 48 |
*/ |
| 49 |
Double_t charge_tracker_rigidity_wolfgang(Double_t *x, Double_t *par){ |
| 50 |
|
| 51 |
Float_t C0[9] = {0 , 1 , 2 , 3 , 4 , 5 , 6 , 8 , 90 }; |
| 52 |
Float_t B0[9] = {0 , 0.867439 , 4.04174 , 9.10558 , 15.3956 , 20.6338 , 25.5426 , 32.1287 , 1000 }; |
| 53 |
Float_t B1[9] = {0 , -0.128483 , -3.18795 , -3.38137 , -9.61749 , -6.44302 , -8.23092 , -7.30003 , 0 }; |
| 54 |
Float_t B2[9] = {0 , 0.764514 , 19.7088 , 38.8391 , 58.1051 , 58.8905 , 55.1826 , 58.9612 , 0 }; |
| 55 |
Float_t B3[9] = {0 , 0.00176465 , -9.23937 , -28.0735 , -50.4548 , -54.3154 , -50.3964 , -82.3765 , 0 }; |
| 56 |
Float_t B4[9] = {0 , -0.016259 , 1.32851 , 6.68077 , 14.3786 , 16.3259 , 15.835 , 54.1588 , 0 }; |
| 57 |
|
| 58 |
|
| 59 |
Float_t yl,yh,m,b,chelp; |
| 60 |
|
| 61 |
// Float_t x1[9],y1[9]; |
| 62 |
// Int_t n1 = 9; |
| 63 |
|
| 64 |
Float_t charge = 1000.; |
| 65 |
Float_t def = 1./((Float_t)x[0]); |
| 66 |
// WM 05-03-2009 calibration valid only for def>0 |
| 67 |
def = fabs(def); |
| 68 |
Float_t dedxtrk = (Float_t) x[1]; dedxtrk = dedxtrk/1.3; |
| 69 |
|
| 70 |
Float_t ym = dedxtrk; |
| 71 |
Float_t xb = def; |
| 72 |
|
| 73 |
for ( Int_t jj=0; jj<8; jj++ ){ |
| 74 |
yl = B0[jj] + B1[jj]*xb + B2[jj]*xb*xb + B3[jj]*xb*xb*xb + B4[jj]*xb*xb*xb*xb; |
| 75 |
yh = B0[jj+1] + B1[jj+1]*xb + B2[jj+1]*xb*xb + B3[jj+1]*xb*xb*xb + B4[jj+1]*xb*xb*xb*xb; |
| 76 |
|
| 77 |
// if ((yl<ym)&&(ym<yh)){ |
| 78 |
//WM 05-03-2009 dEdx=0 => charge=0 |
| 79 |
if ((yl<=ym)&&(ym<yh)){ |
| 80 |
m = (C0[jj+1]*C0[jj+1] - C0[jj]*C0[jj]) / (yh - yl); |
| 81 |
b = (C0[jj]*C0[jj]) - m*yl; |
| 82 |
chelp = m*ym + b; |
| 83 |
charge= TMath::Sqrt(chelp); |
| 84 |
} |
| 85 |
} |
| 86 |
|
| 87 |
return (Double_t)charge ; |
| 88 |
|
| 89 |
} |
| 90 |
//////////////////////////////////////////////////////////// |
| 91 |
/// INITIALIZATIONS & GLOBAL PARAMETERS |
| 92 |
//////////////////////////////////////////////////////////// |
| 93 |
|
| 94 |
TrkNuclei_parameters * TrkNuclei_parameters::_parameters = 0; |
| 95 |
|
| 96 |
|
| 97 |
/** |
| 98 |
* Delete |
| 99 |
*/ |
| 100 |
void TrkNuclei_parameters::Delete(){ |
| 101 |
if(charge_vs_beta)charge_vs_beta->Delete(); |
| 102 |
if(charge_vs_rig)charge_vs_rig->Delete(); |
| 103 |
_parameters=NULL; |
| 104 |
} |
| 105 |
|
| 106 |
/** |
| 107 |
* Set external parameters to default values |
| 108 |
*/ |
| 109 |
void TrkNuclei_parameters::SetDefault(){ |
| 110 |
|
| 111 |
cout <<endl << "---------------------------------------------"<<endl; |
| 112 |
cout << "TrkNuclei --> SET DEFAULT PARAMETERS" << endl; |
| 113 |
cout << "---------------------------------------------"<<endl; |
| 114 |
|
| 115 |
cout << endl << " Charge vs beta "<<endl; |
| 116 |
cout << " >> Wolfgang <<"<<endl; |
| 117 |
charge_vs_beta = new TF2("fch_vs_beta",charge_tracker_beta_wolfgang,0.,1.1,0.,100.); |
| 118 |
|
| 119 |
cout << endl << " Charge vs rigidity "<<endl; |
| 120 |
charge_vs_rig = new TF2("fch_vs_rig",charge_tracker_rigidity_wolfgang,0.,10000.,0.,100.); |
| 121 |
cout << " >> Wolfgang <<"<<endl; |
| 122 |
cout << "---------------------------------------------"<<endl; |
| 123 |
|
| 124 |
} |
| 125 |
/* |
| 126 |
* Set new charge-vs-beta parameterization. |
| 127 |
* Input function should be a 2D function: the first variable is beta, the second one the dedx. |
| 128 |
*/ |
| 129 |
bool TrkNuclei_parameters::SetZ_Beta(TF2* f){ |
| 130 |
if(!f)return false; |
| 131 |
cout << "TrkNuclei --> Setting new charge-vs-beta parameterization"<<endl; |
| 132 |
if(charge_vs_beta)charge_vs_beta->Delete(); |
| 133 |
charge_vs_beta = f; |
| 134 |
return true; |
| 135 |
}; |
| 136 |
/* |
| 137 |
* Set new charge-vs-beta parameterization. |
| 138 |
* Input function should be a 2D function: the first variable is beta, the second one the dedx. |
| 139 |
*/ |
| 140 |
bool TrkNuclei_parameters::SetZ_Beta( const char* formula ){ |
| 141 |
TF2 *f = new TF2("fch_vs_beta",formula,0.,1.1,0.,100.); |
| 142 |
return SetZ_Beta(f); |
| 143 |
}; |
| 144 |
/* |
| 145 |
* Set new charge-vs-rigidity parameterization. |
| 146 |
* Input function should be a 2D function: the first variable is the rigidity (GV), the second one the dedx. |
| 147 |
*/ |
| 148 |
bool TrkNuclei_parameters::SetZ_Rigidity(TF2* f){ |
| 149 |
if(!f)return false; |
| 150 |
cout << "TrkNuclei --> Setting new charge-vs-rigidity parameterization"<<endl; |
| 151 |
if(charge_vs_rig)charge_vs_rig->Delete(); |
| 152 |
charge_vs_rig = f; |
| 153 |
return true; |
| 154 |
}; |
| 155 |
/* |
| 156 |
* Set new charge-vs-rigidity parameterization. |
| 157 |
* Input function should be a 2D function: the first variable is the rigidity (GV), the second one the dedx. |
| 158 |
*/ |
| 159 |
bool TrkNuclei_parameters::SetZ_Rigidity( const char* formula ){ |
| 160 |
TF2 *f = new TF2("fch_vs_rig",formula,0.,1.1,0.,100.); |
| 161 |
return SetZ_Rigidity(f); |
| 162 |
}; |
| 163 |
|
| 164 |
//////////////////////////////////////////////////////////// |
| 165 |
/// CLASS IMPLEMENTATION |
| 166 |
//////////////////////////////////////////////////////////// |
| 167 |
|
| 168 |
/** |
| 169 |
* Set event |
| 170 |
*/ |
| 171 |
bool TrkNuclei::Set(TrkTrack* trk){ |
| 172 |
if(!trk)return false; |
| 173 |
if(!nutrk)nutrk = new TrkTrack(); |
| 174 |
*nutrk = *trk; |
| 175 |
// for(int i=0; i<6; i++){ |
| 176 |
// dedx_x[i]=trk->dedx_x[i]; |
| 177 |
// dedx_y[i]=trk->dedx_y[i]; |
| 178 |
// maxs_x[i]=trk->GetClusterX_MaxStrip(i); |
| 179 |
// maxs_y[i]=trk->GetClusterY_MaxStrip(i); |
| 180 |
// } |
| 181 |
return true; |
| 182 |
}; |
| 183 |
bool TrkNuclei::Set(PamTrack* trk){ |
| 184 |
if(!trk)return false; |
| 185 |
return Set(trk->GetTrkTrack()); |
| 186 |
}; |
| 187 |
bool TrkNuclei::Set(PamLevel2* l2,int ntr){ |
| 188 |
if(!l2)return false; |
| 189 |
if(!l2->GetTrkLevel2())return false; |
| 190 |
if(ntr>l2->GetTrkLevel2()->GetNTracks()||ntr<0)return false; |
| 191 |
return Set(l2->GetTrack(ntr)->GetTrkTrack()); |
| 192 |
}; |
| 193 |
/** |
| 194 |
* Reset event |
| 195 |
*/ |
| 196 |
void TrkNuclei::Reset(){ |
| 197 |
if(nutrk)nutrk->Clear(); |
| 198 |
// for(int i=0; i<6; i++){ |
| 199 |
// dedx_x[i]=0.; |
| 200 |
// dedx_y[i]=0.; |
| 201 |
// maxs_x[i]=-1; |
| 202 |
// maxs_y[i]=-1; |
| 203 |
// } |
| 204 |
} |
| 205 |
/** |
| 206 |
* Method to evaluate the particle charge in each view, as a function of the particle velocity. |
| 207 |
* @param ip plane (0-5) |
| 208 |
* @param iv view (0=x 1=y) |
| 209 |
* @param beta particle velocity |
| 210 |
*/ |
| 211 |
Float_t TrkNuclei::GetZ_Beta(int ip, int iv, float beta){ |
| 212 |
|
| 213 |
if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.; |
| 214 |
return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX(ip,iv)); |
| 215 |
|
| 216 |
} |
| 217 |
/** |
| 218 |
* Method to evaluate the particle charge in each view, as a function of the particle velocity. |
| 219 |
* @param ip plane (0-5) |
| 220 |
* @param beta particle velocity |
| 221 |
*/ |
| 222 |
Float_t TrkNuclei::GetZ_Beta(int ip, float beta){ |
| 223 |
|
| 224 |
if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.; |
| 225 |
return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX(ip)); |
| 226 |
|
| 227 |
} |
| 228 |
/** |
| 229 |
* Method to evaluate the average particle charge, as a function of the particle velocity. |
| 230 |
* @param beta particle velocity |
| 231 |
*/ |
| 232 |
Float_t TrkNuclei::GetZ_Beta(float beta){ |
| 233 |
|
| 234 |
if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.; |
| 235 |
return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX()); |
| 236 |
|
| 237 |
} |
| 238 |
/** |
| 239 |
* Method to evaluate the particle charge in each view, as a function of the particle rigidity. |
| 240 |
* @param ip plane (0-5) |
| 241 |
* @param iv view (0=x 1=y) |
| 242 |
* @param rig particle rigidity |
| 243 |
*/ |
| 244 |
Float_t TrkNuclei::GetZ_Rigidity(int ip, int iv, float rig){ |
| 245 |
|
| 246 |
if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.; |
| 247 |
return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX(ip,iv)); |
| 248 |
|
| 249 |
} |
| 250 |
/** |
| 251 |
* Method to evaluate the particle charge in each view, as a function of the particle rigidity. |
| 252 |
* @param ip plane (0-5) |
| 253 |
* @param rig particle rigidity |
| 254 |
*/ |
| 255 |
Float_t TrkNuclei::GetZ_Rigidity(int ip, float rig){ |
| 256 |
|
| 257 |
if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.; |
| 258 |
return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX(ip)); |
| 259 |
|
| 260 |
} |
| 261 |
/** |
| 262 |
* Method to evaluate the particle charge in each view, as a function of the particle rigidity. |
| 263 |
* @param rig particle rigidity |
| 264 |
*/ |
| 265 |
Float_t TrkNuclei::GetZ_Rigidity(float rig){ |
| 266 |
|
| 267 |
if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.; |
| 268 |
return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX()); |
| 269 |
|
| 270 |
} |
| 271 |
|
| 272 |
ClassImp(TrkNuclei); |
Parent Directory
| 
Revision Log