/[PAMELA software]/tracker/flight/TrkNuclei/src/TrkNuclei.cpp
ViewVC logotype

Contents of /tracker/flight/TrkNuclei/src/TrkNuclei.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Tue Jan 20 10:29:57 2009 UTC (15 years, 10 months ago) by pam-fi
Branch: tracker
CVS Tags: V00
Changes since 1.1: +0 -0 lines
Error occurred while calculating annotation data.
tracker utilities

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 Float_t dedxtrk = (Float_t) x[1]; dedxtrk = dedxtrk/1.3;
67
68 Float_t ym = dedxtrk;
69 Float_t xb = def;
70
71 for ( Int_t jj=0; jj<8; jj++ ){
72 yl = B0[jj] + B1[jj]*xb + B2[jj]*xb*xb + B3[jj]*xb*xb*xb + B4[jj]*xb*xb*xb*xb;
73 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;
74
75 if ((yl<ym)&&(ym<yh)){
76 m = (C0[jj+1]*C0[jj+1] - C0[jj]*C0[jj]) / (yh - yl);
77 b = (C0[jj]*C0[jj]) - m*yl;
78 chelp = m*ym + b;
79 charge= TMath::Sqrt(chelp);
80 }
81 }
82
83 return (Double_t)charge ;
84
85 }
86 ////////////////////////////////////////////////////////////
87 /// INITIALIZATIONS & GLOBAL PARAMETERS
88 ////////////////////////////////////////////////////////////
89
90 TrkNuclei_parameters * TrkNuclei_parameters::_parameters = 0;
91
92
93 /**
94 * Delete
95 */
96 void TrkNuclei_parameters::Delete(){
97 if(charge_vs_beta)charge_vs_beta->Delete();
98 if(charge_vs_rig)charge_vs_rig->Delete();
99 _parameters=NULL;
100 }
101
102 /**
103 * Set external parameters to default values
104 */
105 void TrkNuclei_parameters::SetDefault(){
106
107 cout <<endl << "---------------------------------------------"<<endl;
108 cout << "TrkNuclei --> SET DEFAULT PARAMETERS" << endl;
109 cout << "---------------------------------------------"<<endl;
110
111 cout << endl << " Charge vs beta "<<endl;
112 cout << " >> Wolfgang <<"<<endl;
113 charge_vs_beta = new TF2("fch_vs_beta",charge_tracker_beta_wolfgang,0.,1.1,0.,100.);
114
115 cout << endl << " Charge vs rigidity "<<endl;
116 charge_vs_rig = new TF2("fch_vs_rig",charge_tracker_rigidity_wolfgang,0.,10000.,0.,100.);
117 cout << " >> Wolfgang <<"<<endl;
118 cout << "---------------------------------------------"<<endl;
119
120 }
121 /*
122 * Set new charge-vs-beta parameterization.
123 * Input function should be a 2D function: the first variable is beta, the second one the dedx.
124 */
125 bool TrkNuclei_parameters::SetZ_Beta(TF2* f){
126 if(!f)return false;
127 cout << "TrkNuclei --> Setting new charge-vs-beta parameterization"<<endl;
128 if(charge_vs_beta)charge_vs_beta->Delete();
129 charge_vs_beta = f;
130 return true;
131 };
132 /*
133 * Set new charge-vs-beta parameterization.
134 * Input function should be a 2D function: the first variable is beta, the second one the dedx.
135 */
136 bool TrkNuclei_parameters::SetZ_Beta( const char* formula ){
137 TF2 *f = new TF2("fch_vs_beta",formula,0.,1.1,0.,100.);
138 return SetZ_Beta(f);
139 };
140 /*
141 * Set new charge-vs-rigidity parameterization.
142 * Input function should be a 2D function: the first variable is the rigidity (GV), the second one the dedx.
143 */
144 bool TrkNuclei_parameters::SetZ_Rigidity(TF2* f){
145 if(!f)return false;
146 cout << "TrkNuclei --> Setting new charge-vs-rigidity parameterization"<<endl;
147 if(charge_vs_rig)charge_vs_rig->Delete();
148 charge_vs_rig = f;
149 return true;
150 };
151 /*
152 * Set new charge-vs-rigidity parameterization.
153 * Input function should be a 2D function: the first variable is the rigidity (GV), the second one the dedx.
154 */
155 bool TrkNuclei_parameters::SetZ_Rigidity( const char* formula ){
156 TF2 *f = new TF2("fch_vs_rig",formula,0.,1.1,0.,100.);
157 return SetZ_Rigidity(f);
158 };
159
160 ////////////////////////////////////////////////////////////
161 /// CLASS IMPLEMENTATION
162 ////////////////////////////////////////////////////////////
163
164 /**
165 * Set event
166 */
167 bool TrkNuclei::Set(TrkTrack* trk){
168 if(!trk)return false;
169 if(!nutrk)nutrk = new TrkTrack();
170 *nutrk = *trk;
171 // for(int i=0; i<6; i++){
172 // dedx_x[i]=trk->dedx_x[i];
173 // dedx_y[i]=trk->dedx_y[i];
174 // maxs_x[i]=trk->GetClusterX_MaxStrip(i);
175 // maxs_y[i]=trk->GetClusterY_MaxStrip(i);
176 // }
177 return true;
178 };
179 bool TrkNuclei::Set(PamTrack* trk){
180 if(!trk)return false;
181 return Set(trk->GetTrkTrack());
182 };
183 bool TrkNuclei::Set(PamLevel2* l2,int ntr){
184 if(!l2)return false;
185 if(!l2->GetTrkLevel2())return false;
186 if(ntr>l2->GetTrkLevel2()->GetNTracks()||ntr<0)return false;
187 return Set(l2->GetTrack(ntr)->GetTrkTrack());
188 };
189 /**
190 * Reset event
191 */
192 void TrkNuclei::Reset(){
193 if(nutrk)nutrk->Clear();
194 // for(int i=0; i<6; i++){
195 // dedx_x[i]=0.;
196 // dedx_y[i]=0.;
197 // maxs_x[i]=-1;
198 // maxs_y[i]=-1;
199 // }
200 }
201 /**
202 * Method to evaluate the particle charge in each view, as a function of the particle velocity.
203 * @param ip plane (0-5)
204 * @param iv view (0=x 1=y)
205 * @param beta particle velocity
206 */
207 Float_t TrkNuclei::GetZ_Beta(int ip, int iv, float beta){
208
209 if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.;
210 return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX(ip,iv));
211
212 }
213 /**
214 * Method to evaluate the particle charge in each view, as a function of the particle velocity.
215 * @param ip plane (0-5)
216 * @param beta particle velocity
217 */
218 Float_t TrkNuclei::GetZ_Beta(int ip, float beta){
219
220 if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.;
221 return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX(ip));
222
223 }
224 /**
225 * Method to evaluate the average particle charge, as a function of the particle velocity.
226 * @param beta particle velocity
227 */
228 Float_t TrkNuclei::GetZ_Beta(float beta){
229
230 if( !TrkNuclei_parameters::Get()->GetZ_Beta() )return 0.;
231 return TrkNuclei_parameters::Get()->GetZ_Beta()->Eval(beta,nutrk->GetDEDX());
232
233 }
234 /**
235 * Method to evaluate the particle charge in each view, as a function of the particle rigidity.
236 * @param ip plane (0-5)
237 * @param iv view (0=x 1=y)
238 * @param rig particle rigidity
239 */
240 Float_t TrkNuclei::GetZ_Rigidity(int ip, int iv, float rig){
241
242 if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.;
243 return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX(ip,iv));
244
245 }
246 /**
247 * Method to evaluate the particle charge in each view, as a function of the particle rigidity.
248 * @param ip plane (0-5)
249 * @param rig particle rigidity
250 */
251 Float_t TrkNuclei::GetZ_Rigidity(int ip, float rig){
252
253 if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.;
254 return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX(ip));
255
256 }
257 /**
258 * Method to evaluate the particle charge in each view, as a function of the particle rigidity.
259 * @param rig particle rigidity
260 */
261 Float_t TrkNuclei::GetZ_Rigidity(float rig){
262
263 if( !TrkNuclei_parameters::Get()->GetZ_Rigidity() )return 0.;
264 return TrkNuclei_parameters::Get()->GetZ_Rigidity()->Eval(rig,nutrk->GetDEDX());
265
266 }
267
268 ClassImp(TrkNuclei);

  ViewVC Help
Powered by ViewVC 1.1.23