11 |
extern "C" { |
extern "C" { |
12 |
|
|
13 |
// int readetaparam_(); |
// int readetaparam_(); |
14 |
float cog_(int*,int*); |
float cog_(int*,int*); |
15 |
float pfaeta_(int*,float*); |
float pfaeta_(int*,float*); |
16 |
float pfaeta2_(int*,float*); |
float pfaeta2_(int*,float*); |
17 |
float pfaeta3_(int*,float*); |
float pfaeta3_(int*,float*); |
18 |
float pfaeta4_(int*,float*); |
float pfaeta4_(int*,float*); |
19 |
|
float pfaetal_(int*,float*); |
20 |
|
float digsat_(int*); |
21 |
|
int npfastrips_(int*,float*); |
22 |
|
|
23 |
} |
} |
24 |
//-------------------------------------- |
//-------------------------------------- |
27 |
//-------------------------------------- |
//-------------------------------------- |
28 |
TrkCluster::TrkCluster(){ |
TrkCluster::TrkCluster(){ |
29 |
|
|
30 |
view = 0; |
// cout << "TrkCluster::TrkCluster()"<<endl; |
31 |
maxs = 0; |
view = -1; |
32 |
indmax = 0; |
maxs = -1; |
33 |
|
indmax = -1; |
34 |
CLlength = 0; |
|
35 |
clsignal = 0; |
CLlength = 0; |
36 |
clsigma = 0; |
clsignal = 0; |
37 |
cladc = 0; |
clsigma = 0; |
38 |
clbad = 0; |
cladc = 0; |
39 |
|
clbad = 0; |
40 |
|
|
41 |
}; |
}; |
42 |
//-------------------------------------- |
//-------------------------------------- |
43 |
// |
// |
44 |
// |
// |
45 |
//-------------------------------------- |
//-------------------------------------- |
|
TrkCluster::~TrkCluster(){ |
|
|
|
|
|
if(CLlength){ |
|
|
delete [] clsignal; |
|
|
delete [] clsigma; |
|
|
delete [] cladc; |
|
|
delete [] clbad; |
|
|
} |
|
|
}; |
|
|
//-------------------------------------- |
|
|
// |
|
|
// |
|
|
//-------------------------------------- |
|
46 |
TrkCluster::TrkCluster(const TrkCluster& t){ |
TrkCluster::TrkCluster(const TrkCluster& t){ |
47 |
|
|
48 |
view = t.view; |
view = t.view; |
49 |
maxs = t.maxs; |
maxs = t.maxs; |
50 |
indmax = t.indmax; |
indmax = t.indmax; |
51 |
|
|
52 |
CLlength = t.CLlength; |
CLlength = t.CLlength; |
53 |
if(CLlength){ |
if(CLlength){ |
54 |
clsignal = new Float_t[CLlength]; |
clsignal = new Float_t[CLlength]; |
55 |
clsigma = new Float_t[CLlength]; |
clsigma = new Float_t[CLlength]; |
56 |
cladc = new Int_t[CLlength]; |
cladc = new Int_t[CLlength]; |
57 |
clbad = new Bool_t[CLlength]; |
clbad = new Bool_t[CLlength]; |
58 |
for(Int_t i=0; i<CLlength;i++){ |
for(Int_t i=0; i<CLlength;i++){ |
59 |
clsignal[i] = t.clsignal[i]; |
clsignal[i] = t.clsignal[i]; |
60 |
clsigma[i] = t.clsigma[i]; |
clsigma[i] = t.clsigma[i]; |
61 |
cladc[i] = t.cladc[i]; |
cladc[i] = t.cladc[i]; |
62 |
clbad[i] = t.clbad[i]; |
clbad[i] = t.clbad[i]; |
|
}; |
|
63 |
}; |
}; |
64 |
|
}; |
65 |
}; |
}; |
66 |
//-------------------------------------- |
//-------------------------------------- |
67 |
// |
// |
68 |
// |
// |
69 |
//-------------------------------------- |
//-------------------------------------- |
70 |
void TrkCluster::Clear(){ |
void TrkCluster::Clear(){ |
71 |
|
|
72 |
view = 0; |
// cout << "void TrkCluster::Clear()"<<endl; |
73 |
maxs = 0; |
if(CLlength){ |
74 |
indmax = 0; |
delete [] clsignal; |
75 |
|
delete [] clsigma; |
76 |
CLlength = 0; |
delete [] cladc; |
77 |
clsignal = 0; |
delete [] clbad; |
78 |
clsigma = 0; |
} |
79 |
cladc = 0; |
|
80 |
clbad = 0; |
view = 0; |
81 |
|
maxs = 0; |
82 |
|
indmax = 0; |
83 |
|
|
84 |
|
CLlength = 0; |
85 |
|
clsignal = 0; |
86 |
|
clsigma = 0; |
87 |
|
cladc = 0; |
88 |
|
clbad = 0; |
89 |
|
|
90 |
}; |
}; |
91 |
//-------------------------------------- |
//-------------------------------------- |
97 |
* strips, around maxs, having a significant signal. |
* strips, around maxs, having a significant signal. |
98 |
* @param nstrip Maximum number of strips. |
* @param nstrip Maximum number of strips. |
99 |
* @param cut Inclusion cut ( s > cut*sigma ). |
* @param cut Inclusion cut ( s > cut*sigma ). |
100 |
|
* @param force Falg to force the PFA strip-inclusion pattern (nstrip>0) |
101 |
* If nstrip<=0 only the inclusion cut is used to determine the cluster size. |
* If nstrip<=0 only the inclusion cut is used to determine the cluster size. |
102 |
*/ |
*/ |
103 |
Float_t TrkCluster::GetSignal(Int_t nstrip, Float_t cut){ |
Float_t TrkCluster::GetSignal(Int_t nstrip, Float_t cut, Bool_t force){ |
104 |
|
|
105 |
|
if(CLlength<=0)return 0; |
106 |
|
|
107 |
Float_t s = 0; |
Float_t s = 0; |
108 |
|
|
109 |
|
//----------------------------------- |
110 |
|
// inlcude strips with s > cut*sigma |
111 |
|
//----------------------------------- |
112 |
|
|
113 |
if( nstrip<=0 ){ |
if( nstrip<=0 ){ |
114 |
// for(Int_t is = 0; is < CLlength; is++){ |
// for(Int_t is = 0; is < CLlength; is++){ |
115 |
// Float_t scut = cut*clsigma[is]; |
// Float_t scut = cut*clsigma[is]; |
128 |
return s; |
return s; |
129 |
}; |
}; |
130 |
|
|
131 |
|
//--------------------------------------------------- |
132 |
|
// inlcude strips with s > cut*sigma, up to nstrip. |
133 |
|
// strips are included in order of decreasing signal |
134 |
|
//--------------------------------------------------- |
135 |
|
if( !force ){ |
136 |
|
|
137 |
|
Int_t il = indmax; |
138 |
|
Int_t ir = indmax; |
139 |
|
Int_t inc = 0; |
140 |
|
|
141 |
|
if( clsignal[indmax] < cut*clsigma[indmax] ) return 0; |
142 |
|
|
143 |
|
while ( inc < nstrip ){ |
144 |
|
Float_t sl = -100000; |
145 |
|
Float_t sr = -100000; |
146 |
|
if( il >= 0 ) sl = clsignal[il]; |
147 |
|
if( ir < CLlength ) sr = clsignal[ir]; |
148 |
|
if( sl == sr && inc == 0 ){ |
149 |
|
s += clsignal[il]; //cout << inc<<" - "<< clsignal[il]<<" "<<s<<endl; |
150 |
|
il--; |
151 |
|
ir++; |
152 |
|
}else if ( sl >= sr && sl > cut*clsigma[il] && inc !=0 ){ |
153 |
|
s += sl;//cout << inc<<" - "<< clsignal[il]<<" "<<s<<endl; |
154 |
|
il--; |
155 |
|
}else if ( sl < sr && sr > cut*clsigma[ir] ){ |
156 |
|
s += sr;//cout << inc<<" - " << clsignal[ir]<<" "<<s<<endl; |
157 |
|
ir++; |
158 |
|
}else break; |
159 |
|
|
160 |
|
inc++; |
161 |
|
} |
162 |
|
return s; |
163 |
|
|
164 |
|
}else{ |
165 |
|
//--------------------------------------------------- |
166 |
|
// evaluate signal using a fixed number of strips, |
167 |
|
// following the PFA inclusion patters |
168 |
|
//--------------------------------------------------- |
169 |
|
// --> signal of the central strip |
170 |
|
Float_t sc = clsignal[indmax]; |
171 |
|
// signal of adjacent strips |
172 |
|
Float_t sl1 = -9999.; |
173 |
|
Float_t sl2 = -9999.; |
174 |
|
Float_t sr1 = -9999.; |
175 |
|
Float_t sr2 = -9999.; |
176 |
|
if(indmax-1>=0) sl1 = clsignal[indmax-1]; |
177 |
|
if(indmax-2>=0) sl2 = clsignal[indmax-2]; |
178 |
|
if(indmax+1<CLlength) sr1 = clsignal[indmax+1]; |
179 |
|
if(indmax+2<CLlength) sr2 = clsignal[indmax+2]; |
180 |
|
|
181 |
|
if(nstrip==1){ |
182 |
|
s = sc; |
183 |
|
}else if(nstrip==2){ |
184 |
|
if( sl1>sr1 && sl1+sc!=0 )s = (sl1+sc); |
185 |
|
if( sl1<sr1 && sr1+sc!=0 )s = (sc+sr1); |
186 |
|
if( sl1==sr1 && sl1 != -9999.){ |
187 |
|
if( clsigma[indmax-1] < clsigma[indmax+1] && sl1+sc!=0 )s = (sl1+sc); |
188 |
|
if( clsigma[indmax-1] > clsigma[indmax+1] && sc+sr1!=0 )s = (sc+sr1); |
189 |
|
} |
190 |
|
}else if(nstrip==3){ |
191 |
|
s = (sl1+sc+sr1); |
192 |
|
}else if(nstrip==4){ |
193 |
|
if( sl2>sr2 && sl2+sl1+sc+sr1!=0 )s = (sl2+sl1+sc+sr1); |
194 |
|
if( sl2<sr2 && sl1+sc+sr1+sr2!=0 )s = (sl1+sc+sr1+sr2); |
195 |
|
if( sl2==sr2 && sl2 != -9999.){ |
196 |
|
if( clsigma[indmax-2] < clsigma[indmax+2] && sl2+sl1+sc+sr1!=0 )s = (sl2+sl1+sc+sr1); |
197 |
|
if( clsigma[indmax-2] > clsigma[indmax+2] && sl1+sc+sr1+sr2!=0 )s = (sl1+sc+sr1+sr2); |
198 |
|
} |
199 |
|
}else if(nstrip==5){ |
200 |
|
s = (sl1+sc+sr1); |
201 |
|
if(sl2 != -9999.)s += sl2; |
202 |
|
if(sr2 != -9999.)s += sr2; |
203 |
|
}else{ |
204 |
|
cout << "Float_t TrkCluster::GetSignal("<<nstrip<<","<<cut<<","<<force<<")- not implemented"<<endl; |
205 |
|
} |
206 |
|
|
207 |
Int_t il = indmax; |
} |
|
Int_t ir = indmax; |
|
|
Int_t inc = 0; |
|
208 |
|
|
209 |
if( clsignal[indmax] < cut*clsigma[indmax] ) return 0; |
return 0.; |
210 |
|
|
|
while ( inc < nstrip ){ |
|
|
Float_t sl = -100000; |
|
|
Float_t sr = -100000; |
|
|
if( il >= 0 ) sl = clsignal[il]; |
|
|
if( ir < CLlength ) sr = clsignal[ir]; |
|
|
if( sl == sr && inc == 0 ){ |
|
|
s += clsignal[il]; //cout << inc<<" - "<< clsignal[il]<<" "<<s<<endl; |
|
|
il--; |
|
|
ir++; |
|
|
}else if ( sl >= sr && sl > cut*clsigma[il] && inc !=0 ){ |
|
|
s += sl;//cout << inc<<" - "<< clsignal[il]<<" "<<s<<endl; |
|
|
il--; |
|
|
}else if ( sl < sr && sr > cut*clsigma[ir] ){ |
|
|
s += sr;//cout << inc<<" - " << clsignal[ir]<<" "<<s<<endl; |
|
|
ir++; |
|
|
}else break; |
|
|
|
|
|
inc++; |
|
|
} |
|
|
return s; |
|
211 |
}; |
}; |
212 |
|
|
213 |
|
|
214 |
/** |
/** |
215 |
including a ( maximum ) fixed number of adjacent strips (with s>0) around the maxs. |
* Evaluate the cluster signal-to-noise, as defined by Turchetta, including a |
216 |
* @param nstrip Number of strips. |
* maximum number of adjacent strips, around maxs, having a significant signal. |
|
*/ |
|
|
/** |
|
|
* Evaluate the cluster signal-to-noise, as defined by Turchetta, including a maximum number of adjacent strips, around maxs, having a significant signal. |
|
217 |
* @param nstrip Maximum number of strips. |
* @param nstrip Maximum number of strips. |
218 |
* @param cut Inclusion cut ( s > cut*sigma ). |
* @param cut Inclusion cut ( s > cut*sigma ). |
219 |
* If nstrip<=0 only the inclusion cut is used to determine the cluster size. |
* If nstrip<=0 only the inclusion cut is used to determine the cluster size. |
220 |
*/ |
*/ |
221 |
Float_t TrkCluster::GetSignalToNoise(Int_t nstrip, Float_t cut){ |
Float_t TrkCluster::GetSignalToNoise(Int_t nstrip, Float_t cut){ |
222 |
|
|
223 |
|
if(CLlength<=0)return 0; |
224 |
|
|
225 |
Float_t sn = 0; |
Float_t sn = 0; |
226 |
|
|
227 |
if( nstrip<=0 ){ |
if( nstrip<=0 ){ |
271 |
* @param cut Inclusion cut. |
* @param cut Inclusion cut. |
272 |
*/ |
*/ |
273 |
Int_t TrkCluster::GetMultiplicity(Float_t cut){ |
Int_t TrkCluster::GetMultiplicity(Float_t cut){ |
274 |
Int_t m = 0; |
|
275 |
for(Int_t is = 0; is < CLlength; is++){ |
if(CLlength<=0)return 0; |
276 |
Float_t scut = cut*clsigma[is]; |
|
277 |
if(clsignal[is] > scut) m++; |
Int_t m = 0; |
278 |
}; |
|
279 |
return m; |
for(Int_t is = indmax+1; is < CLlength; is++){ |
280 |
|
Float_t scut = cut*clsigma[is]; |
281 |
|
if(clsignal[is] > scut) m++; |
282 |
|
else break; |
283 |
|
}; |
284 |
|
for(Int_t is = indmax; is >=0; is--){ |
285 |
|
Float_t scut = cut*clsigma[is]; |
286 |
|
if(clsignal[is] > scut) m++; |
287 |
|
else break; |
288 |
|
}; |
289 |
|
return m; |
290 |
}; |
}; |
291 |
/** |
/** |
292 |
* True if the cluster contains bad strips. |
* True if the cluster contains bad strips. |
293 |
* @param nbad Number of strips around the maximum. |
* @param nbad Number of strips around the maximum. |
294 |
*/ |
*/ |
295 |
Bool_t TrkCluster::IsBad(Int_t nbad){ |
Bool_t TrkCluster::IsBad(Int_t nbad){ |
296 |
|
|
297 |
/* Float_t max = 0; |
if(CLlength<=0)return 0; |
298 |
Int_t imax = 0; |
|
299 |
for(Int_t is = 0; is < CLlength; is++){ |
Int_t il,ir; |
300 |
if(clsignal[is] > max){ |
il = indmax; |
301 |
max = clsignal[is]; |
ir = indmax; |
302 |
imax = is; |
for(Int_t i=1; i<nbad; i++){ |
303 |
}; |
if (ir == CLlength-1 && il == 0)break; |
304 |
}; |
else if (ir == CLlength-1 && il != 0)il--; |
305 |
|
else if (ir != CLlength-1 && il == 0)ir++; |
306 |
Int_t il,ir; |
else{ |
307 |
il = imax; |
if(clsignal[il-1] > clsignal[ir+1])il--; |
308 |
ir = imax;*/ |
else ir++; |
|
|
|
|
Int_t il,ir; |
|
|
il = indmax; |
|
|
ir = indmax; |
|
|
for(Int_t i=1; i<nbad; i++){ |
|
|
if (ir == CLlength && il == 0)break; |
|
|
else if (ir == CLlength && il != 0)il--; |
|
|
else if (ir != CLlength && il == 0)ir++; |
|
|
else{ |
|
|
if(clsignal[il-1] > clsignal[ir+1])il--; |
|
|
else ir++; |
|
|
} |
|
309 |
} |
} |
310 |
Int_t isbad = 0; |
} |
311 |
for(Int_t i=il; i<=ir; i++)isbad += clbad[i]; |
Int_t isbad = 0; |
312 |
|
for(Int_t i=il; i<=ir; i++)isbad += clbad[i]; |
313 |
return ( isbad != nbad ); |
|
314 |
|
return ( isbad != nbad ); |
315 |
}; |
}; |
316 |
|
/** |
317 |
|
* True if the cluster contains saturated strips. |
318 |
|
* @param nbad Number of strips around the maximum. |
319 |
|
*/ |
320 |
|
Bool_t TrkCluster::IsSaturated(Int_t nbad){ |
321 |
|
|
322 |
|
if(CLlength<=0)return 0; |
323 |
|
|
324 |
|
Int_t il,ir; |
325 |
|
il = indmax; |
326 |
|
ir = indmax; |
327 |
|
for(Int_t i=1; i<nbad; i++){ |
328 |
|
if (ir == CLlength-1 && il == 0)break; |
329 |
|
else if (ir == CLlength-1 && il != 0)il--; |
330 |
|
else if (ir != CLlength-1 && il == 0)ir++; |
331 |
|
else{ |
332 |
|
if(clsignal[il-1] > clsignal[ir+1])il--; |
333 |
|
else ir++; |
334 |
|
} |
335 |
|
} |
336 |
|
Int_t isbad = 0; |
337 |
|
for(Int_t i=il; i<=ir; i++){ |
338 |
|
if( IsX() && cladc[i] > 2980 )isbad++; |
339 |
|
if( IsY() && cladc[i] < 80 )isbad++; |
340 |
|
} |
341 |
|
return ( isbad != 0 ); |
342 |
|
|
343 |
|
} |
344 |
//-------------------------------------- |
//-------------------------------------- |
345 |
// |
// |
346 |
// |
// |
347 |
//-------------------------------------- |
//-------------------------------------- |
348 |
void TrkCluster::Dump(){ |
void TrkCluster::Dump(){ |
349 |
|
|
350 |
cout << "----- Cluster" << endl; |
cout << "----- Cluster" << endl; |
351 |
cout << "View "<<view << " - Ladder "<<GetLadder()<<endl; |
cout << "View "<<view << " - Ladder "<<GetLadder()<<endl; |
352 |
cout << "Position of maximun "<< maxs <<endl; |
cout << "Position of maximun "<< maxs <<endl; |
353 |
cout << "Multiplicity "<< GetMultiplicity() <<endl; |
cout << "Multiplicity "<< GetMultiplicity() <<endl; |
354 |
cout << "Tot signal "<< GetSignal() << " (ADC channels)"<<endl ; |
cout << "Tot signal "<< GetSignal() << " (ADC channels)"<<endl ; |
355 |
cout << "Signal/Noise "<< GetSignalToNoise(); |
cout << "Signal/Noise "<< GetSignalToNoise()<<endl; |
356 |
cout <<endl<< "Strip signals "; |
cout << "COG "<< GetCOG(0)<<endl;; |
357 |
for(Int_t i =0; i<CLlength; i++)cout << " " <<clsignal[i]; |
cout << "Strip signals "; |
358 |
cout <<endl<< "Strip sigmas "; |
for(Int_t i =0; i<CLlength; i++)cout << " " <<clsignal[i]; |
359 |
for(Int_t i =0; i<CLlength; i++)cout << " " <<clsigma[i]; |
cout <<endl<< "Strip sigmas "; |
360 |
cout <<endl<< "Strip ADC "; |
for(Int_t i =0; i<CLlength; i++)cout << " " <<clsigma[i]; |
361 |
for(Int_t i =0; i<CLlength; i++)cout << " " <<cladc[i]; |
cout <<endl<< "Strip ADC "; |
362 |
cout <<endl<< "Strip BAD "; |
for(Int_t i =0; i<CLlength; i++)cout << " " <<cladc[i]; |
363 |
for(Int_t i =0; i<CLlength; i++){ |
cout <<endl<< "Strip BAD "; |
364 |
if(i==indmax)cout << " *" <<clbad[i]<<"*"; |
for(Int_t i =0; i<CLlength; i++){ |
365 |
else cout << " " <<clbad[i]; |
if(i==indmax)cout << " *" <<clbad[i]<<"*"; |
366 |
} |
else cout << " " <<clbad[i]; |
367 |
cout << endl; |
} |
368 |
|
cout << endl; |
369 |
|
|
370 |
} |
} |
371 |
//-------------------------------------- |
//-------------------------------------- |
375 |
/** |
/** |
376 |
* Method to fill a level1 struct with only one cluster (done to use F77 p.f.a. routines on a cluster basis). |
* Method to fill a level1 struct with only one cluster (done to use F77 p.f.a. routines on a cluster basis). |
377 |
*/ |
*/ |
378 |
cTrkLevel1* TrkCluster::GetLevel1Struct(){ |
void TrkCluster::GetLevel1Struct(cTrkLevel1* l1){ |
379 |
|
|
380 |
// cTrkLevel1* l1 = new cTrkLevel1; |
// cTrkLevel1* l1 = new cTrkLevel1; |
381 |
|
|
382 |
cTrkLevel1* l1 = &level1event_ ; |
// cTrkLevel1* l1 = &level1event_ ; |
383 |
|
|
384 |
l1->nclstr1 = 1; |
l1->nclstr1 = 1; |
385 |
l1->view[0] = view; |
l1->view[0] = view; |
397 |
l1->clbad[i] = clbad[i]; |
l1->clbad[i] = clbad[i]; |
398 |
}; |
}; |
399 |
|
|
400 |
return l1; |
// return l1; |
401 |
}; |
}; |
402 |
//-------------------------------------- |
//-------------------------------------- |
403 |
// |
// |
449 |
// |
// |
450 |
//-------------------------------------- |
//-------------------------------------- |
451 |
/** |
/** |
452 |
* Evaluates the cluster position, in strips, relative to the strip with the maximum signal (TrkCluster::maxs), by applying the non-linear ETA-algorythm. |
* Evaluates the cluster position, in pitch units, relative to the strip |
453 |
|
* with the maximum signal (TrkCluster::maxs), by applying the non-linear |
454 |
|
* ETA-algorythm. |
455 |
* @param neta Number of strips to evaluate ETA. |
* @param neta Number of strips to evaluate ETA. |
456 |
* @param angle Projected angle between particle track and detector plane. |
* @param angle Projected (effective) angle between particle track and detector plane. |
457 |
|
* @landi flag to apply Landi correction |
458 |
* Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle. |
* Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle. |
459 |
*/ |
*/ |
460 |
Float_t TrkCluster::GetETA(Int_t neta, float angle){ |
Float_t TrkCluster::GetETA(Int_t neta, float angle, bool landi){ |
461 |
|
|
462 |
// cout << "GetETA(neta,angle) "<< neta << " "<< angle; |
// cout << "GetETA(neta,angle) "<< neta << " "<< angle; |
463 |
// LoadPfaParam(); |
// LoadPfaParam(); |
464 |
|
|
465 |
|
TrkParams::Load(4); |
466 |
|
if( !TrkParams::IsLoaded(4) ){ |
467 |
|
cout << "Float_t TrkCluster::GetETA(Int_t neta, float angle, bool landi) --- ERROR --- p.f.a. parameters not loaded"<<endl; |
468 |
|
return 0; |
469 |
|
} |
470 |
|
|
471 |
float ax = angle; |
float ax = angle; |
472 |
int ic = 1; |
int ic = 1; |
473 |
GetLevel1Struct(); |
GetLevel1Struct(); |
474 |
if(neta == 0) return pfaeta_(&ic,&ax); |
if( neta == 0 && !landi) return pfaeta_(&ic,&ax); |
475 |
else if(neta == 2) return pfaeta2_(&ic,&ax); |
else if(neta == 0 && landi ) return pfaetal_(&ic,&ax); |
476 |
else if(neta == 3) return pfaeta3_(&ic,&ax); |
else if(neta == 2 ) return pfaeta2_(&ic,&ax); |
477 |
else if(neta == 4) return pfaeta4_(&ic,&ax); |
else if(neta == 3 ) return pfaeta3_(&ic,&ax); |
478 |
else cout << "ETA"<<neta<<" not implemented\n"; |
else if(neta == 4 ) return pfaeta4_(&ic,&ax); |
479 |
|
else cout << "TrkCluster::GetETA("<<neta<<","<<angle<<","<<landi<<") not implemented\n"; |
480 |
return 0; |
return 0; |
481 |
|
|
482 |
}; |
}; |
483 |
|
|
484 |
|
/** |
485 |
|
* Evaluates the cluster position, in pitch units, relative to the strip |
486 |
|
* with the maximum signal (TrkCluster::maxs), by applying the digital |
487 |
|
* algorithm for saturated clusters. |
488 |
|
* |
489 |
|
* @return The cluster position (0 also if if no saturated strip is found). |
490 |
|
*/ |
491 |
|
Float_t TrkCluster::GetDigSat() { |
492 |
|
|
493 |
|
GetLevel1Struct(); |
494 |
|
int ic = 1; |
495 |
|
return digsat_(&ic); |
496 |
|
|
497 |
|
} |
498 |
|
|
499 |
|
/** |
500 |
|
* Evaluates the cluster position, in pitch unit, relative to the strip with |
501 |
|
* the maximum signal (TrkCluster::maxs), by applying the PFA set as default (see TrkParams). |
502 |
|
* @param angle Projected (effective) angle between particle track and detector plane. |
503 |
|
*/ |
504 |
|
Float_t TrkCluster::GetPositionPU(float angle){ |
505 |
|
|
506 |
|
if ( TrkParams::GetPFA() == 0 )return GetETA(0,angle,false); |
507 |
|
else if( TrkParams::GetPFA() == 2 )return GetETA(2,angle,false); |
508 |
|
else if( TrkParams::GetPFA() == 3 )return GetETA(3,angle,false); |
509 |
|
else if( TrkParams::GetPFA() == 4 )return GetETA(4,angle,false); |
510 |
|
else if( TrkParams::GetPFA() == 5 )return GetETA(0,angle,true); |
511 |
|
else if( TrkParams::GetPFA() == 10 )return GetCOG(0); |
512 |
|
else if( TrkParams::GetPFA() == 11 )return GetCOG(1); |
513 |
|
else if( TrkParams::GetPFA() == 12 )return GetCOG(2); |
514 |
|
else if( TrkParams::GetPFA() == 13 )return GetCOG(3); |
515 |
|
else if( TrkParams::GetPFA() == 14 )return GetCOG(4); |
516 |
|
else cout << " TrkCluster::GetPositionPU(float "<<angle<<") -- WARNING -- PFA="<<TrkParams::GetPFA()<<" not implemented"<<endl; |
517 |
|
|
518 |
|
return 0.; |
519 |
|
|
520 |
|
} |
521 |
|
|
522 |
|
/** |
523 |
|
* Give the number of strip used to evaluate the cluster coordinate |
524 |
|
* according to the p.f.a. |
525 |
|
* It returns 0 when the COG is used (in this case the number of strip used |
526 |
|
* equals the multiplicity). |
527 |
|
*/ |
528 |
|
Int_t TrkCluster::GetPFAstrips(float angle){ |
529 |
|
|
530 |
|
float ax = angle; |
531 |
|
int ic = 1; |
532 |
|
GetLevel1Struct(); |
533 |
|
return npfastrips_(&ic,&ax); |
534 |
|
|
535 |
|
} |
536 |
|
|
537 |
//-------------------------------------- |
//-------------------------------------- |
538 |
// |
// |
539 |
// |
// |
540 |
//-------------------------------------- |
//-------------------------------------- |
541 |
TrkLevel1::TrkLevel1(){ |
TrkLevel1::TrkLevel1(){ |
542 |
|
|
543 |
Cluster = new TClonesArray("TrkCluster"); |
// cout << "TrkLevel1::TrkLevel1()"<<endl; |
544 |
|
// Cluster = new TClonesArray("TrkCluster"); |
545 |
|
Cluster = 0; |
546 |
for(Int_t i=0; i<12 ; i++){ |
for(Int_t i=0; i<12 ; i++){ |
547 |
good[i] = -1; |
good[i] = -1; |
548 |
for(Int_t j=0; j<24 ; j++){ |
for(Int_t j=0; j<24 ; j++){ |
549 |
cn[j][i]=0; |
cn[j][i]=0; |
|
// cnrms[j][i]=0; |
|
550 |
cnn[j][i]=0; |
cnn[j][i]=0; |
551 |
}; |
}; |
552 |
}; |
}; |
553 |
|
// TrkParams::SetTrackingMode(); |
554 |
|
// TrkParams::SetPrecisionFactor(); |
555 |
|
// TrkParams::SetStepMin(); |
556 |
|
TrkParams::SetMiniDefault(); |
557 |
|
TrkParams::SetPFA(); |
558 |
|
} |
559 |
|
//-------------------------------------- |
560 |
|
// |
561 |
|
// |
562 |
|
//-------------------------------------- |
563 |
|
void TrkLevel1::Set(){ |
564 |
|
if(!Cluster)Cluster = new TClonesArray("TrkCluster"); |
565 |
} |
} |
566 |
//-------------------------------------- |
//-------------------------------------- |
567 |
// |
// |
582 |
cout <<endl; |
cout <<endl; |
583 |
} |
} |
584 |
|
|
585 |
|
if(!Cluster)return; |
586 |
TClonesArray &t = *Cluster; |
TClonesArray &t = *Cluster; |
587 |
for(int i=0; i<this->nclstr(); i++) ((TrkCluster *)t[i])->Dump(); |
for(int i=0; i<this->nclstr(); i++) ((TrkCluster *)t[i])->Dump(); |
588 |
|
|
589 |
} |
} |
590 |
|
/** |
591 |
|
* \brief Dump processing status |
592 |
|
*/ |
593 |
|
void TrkLevel1::StatusDump(int view){ |
594 |
|
cout << "DSP n. "<<view+1<<" (level1-)status: "<<hex<<showbase<<good[view]<<dec<<endl; |
595 |
|
}; |
596 |
|
/** |
597 |
|
* \brief Check event status |
598 |
|
* |
599 |
|
* Check the event status, according to a flag-mask given as input. |
600 |
|
* Return true if the view passes the check. |
601 |
|
* |
602 |
|
* @param view View number (0-11) |
603 |
|
* @param flagmask Mask of flags to check (eg. flagmask=0x111 no missing packet, |
604 |
|
* no crc error, no software alarm) |
605 |
|
* |
606 |
|
* @see TrkLevel2 class definition to know how the status flag is defined |
607 |
|
* |
608 |
|
*/ |
609 |
|
Bool_t TrkLevel1::StatusCheck(int view, int flagmask){ |
610 |
|
|
611 |
|
if( view<0 || view >= 12)return false; |
612 |
|
return !(good[view]&flagmask); |
613 |
|
|
614 |
|
}; |
615 |
|
|
616 |
|
|
617 |
//-------------------------------------- |
//-------------------------------------- |
618 |
// |
// |
619 |
// |
// |
623 |
*/ |
*/ |
624 |
void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1, Bool_t full){ |
void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1, Bool_t full){ |
625 |
|
|
626 |
|
// cout << "void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1, Bool_t full)"<<endl; |
627 |
|
|
628 |
|
Clear(); |
629 |
// --------------- |
// --------------- |
630 |
// *** CLUSTER *** |
// *** CLUSTER *** |
631 |
// --------------- |
// --------------- |
632 |
TrkCluster* t_cl = new TrkCluster(); |
TrkCluster* t_cl = new TrkCluster(); |
633 |
|
if(!Cluster)Cluster = new TClonesArray("TrkCluster"); |
634 |
TClonesArray &t = *Cluster; |
TClonesArray &t = *Cluster; |
635 |
for(int i=0; i<l1->nclstr1; i++){ |
for(int i=0; i<l1->nclstr1; i++){ |
636 |
|
|
637 |
t_cl->Clear(); |
t_cl->Clear(); |
638 |
if( full || (!full && l1->whichtrack[i]) ){ |
// if( full || (!full && l1->whichtrack[i]) ){ |
639 |
|
|
640 |
t_cl->view = l1->view[i]; |
t_cl->view = l1->view[i]; |
641 |
t_cl->maxs = l1->maxs[i]; |
t_cl->maxs = l1->maxs[i]; |
642 |
t_cl->indmax = l1->indmax[i] - l1->indstart[i]; |
|
643 |
|
if( full || (!full && l1->whichtrack[i]) ){ |
644 |
|
t_cl->indmax = l1->indmax[i] - l1->indstart[i]; |
645 |
Int_t from = l1->indstart[i] -1; |
Int_t from = l1->indstart[i] -1; |
646 |
Int_t to = l1->totCLlength ; |
Int_t to = l1->totCLlength ; |
647 |
if(i != l1->nclstr1-1)to = l1->indstart[i+1] -1 ; |
if(i != l1->nclstr1-1)to = l1->indstart[i+1] -1 ; |
651 |
t_cl->clsigma = new Float_t[t_cl->CLlength]; |
t_cl->clsigma = new Float_t[t_cl->CLlength]; |
652 |
t_cl->cladc = new Int_t[t_cl->CLlength]; |
t_cl->cladc = new Int_t[t_cl->CLlength]; |
653 |
t_cl->clbad = new Bool_t[t_cl->CLlength]; |
t_cl->clbad = new Bool_t[t_cl->CLlength]; |
654 |
|
|
655 |
Int_t index = 0; |
Int_t index = 0; |
656 |
for(Int_t is = from; is < to; is++ ){ |
for(Int_t is = from; is < to; is++ ){ |
657 |
t_cl->clsignal[index] = (Float_t) l1->clsignal[is]; |
t_cl->clsignal[index] = (Float_t) l1->clsignal[is]; |
683 |
* Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common). |
* Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common). |
684 |
*/ |
*/ |
685 |
|
|
686 |
cTrkLevel1* TrkLevel1::GetLevel1Struct() { |
void TrkLevel1::GetLevel1Struct(cTrkLevel1* l1) { |
687 |
|
|
688 |
cTrkLevel1 *l1=0; |
// cTrkLevel1* l1 = &level1event_ ; |
689 |
// |
|
690 |
for(Int_t i=0; i<12 ; i++){ |
for(Int_t i=0; i<12 ; i++){ |
691 |
l1->good[i] = good[i]; |
l1->good[i] = good[i]; |
692 |
for(Int_t j=0; j<24 ; j++){ |
for(Int_t j=0; j<24 ; j++){ |
693 |
l1->cnev[j][i] = cn[j][i]; |
l1->cnev[j][i] = cn[j][i] ; |
694 |
// l1->cnrmsev[j][i] = cnrms[j][i]; |
l1->cnnev[j][i] = cnn[j][i] ; |
695 |
l1->cnnev[j][i] = cnn[j][i]; |
l1->cnrmsev[j][i] = 0. ; |
696 |
}; |
}; |
697 |
|
l1->fshower[i] = 0; |
698 |
}; |
}; |
699 |
|
|
700 |
// *** CLUSTERS *** |
l1->nclstr1=0; |
701 |
l1->nclstr1 = Cluster->GetEntries(); |
l1->totCLlength=0; |
702 |
for(Int_t i=0;i<l1->nclstr1;i++){ |
Int_t index=0; |
703 |
|
if(Cluster){ |
704 |
l1->view[i] = ((TrkCluster *)Cluster->At(i))->view; |
Int_t i=0; |
705 |
l1->maxs[i] = ((TrkCluster *)Cluster->At(i))->maxs; |
for(Int_t ii=0;ii<Cluster->GetEntries();ii++){ |
706 |
// COMPLETARE // |
TrkCluster *clu = GetCluster(ii); |
707 |
// COMPLETARE // |
// ---------------------------------------- |
708 |
// COMPLETARE // |
// attenzione!! |
709 |
// COMPLETARE // |
// se il cluster non e` salvato (view = 0) |
710 |
// COMPLETARE // |
// DEVE essere escluso dal common F77 |
711 |
// COMPLETARE // |
// ---------------------------------------- |
712 |
|
if(clu->view != 0 ){ |
713 |
} |
l1->view[i] = clu->view; |
714 |
// COMPLETARE // |
l1->ladder[i] = clu->GetLadder(); |
715 |
// COMPLETARE // |
l1->maxs[i] = clu->maxs; |
716 |
// COMPLETARE // |
l1->mult[i] = clu->GetMultiplicity(); |
717 |
// COMPLETARE // |
l1->dedx[i] = clu->GetSignal(); |
718 |
// COMPLETARE // |
l1->indstart[i] = index+1; |
719 |
// COMPLETARE // |
l1->indmax[i] = l1->indstart[i] + clu->indmax; |
720 |
return l1; |
l1->totCLlength += clu->CLlength; |
721 |
|
for(Int_t iw=0; iw < clu->CLlength; iw++){ |
722 |
|
l1->clsignal[index] = clu->clsignal[iw]; |
723 |
|
l1->clsigma[index] = clu->clsigma[iw]; |
724 |
|
l1->cladc[index] = clu->cladc[iw]; |
725 |
|
l1->clbad[index] = clu->clbad[iw]; |
726 |
|
index++; |
727 |
|
} |
728 |
|
i++; |
729 |
|
} |
730 |
|
} |
731 |
|
l1->nclstr1 = i; |
732 |
|
} |
733 |
|
|
734 |
|
// return l1; |
735 |
} |
} |
736 |
//-------------------------------------- |
//-------------------------------------- |
737 |
// |
// |
743 |
good[i] = -1; |
good[i] = -1; |
744 |
for(Int_t j=0; j<24 ; j++){ |
for(Int_t j=0; j<24 ; j++){ |
745 |
cn[j][i] = 0; |
cn[j][i] = 0; |
|
// cnrms[j][i] = 0; |
|
746 |
cnn[j][i] = 0; |
cnn[j][i] = 0; |
747 |
}; |
}; |
748 |
}; |
}; |
749 |
// |
// if(Cluster)Cluster->Clear("C"); |
750 |
Cluster->Clear(); |
if(Cluster)Cluster->Delete(); |
751 |
|
|
752 |
} |
} |
753 |
//-------------------------------------- |
//-------------------------------------- |
755 |
// |
// |
756 |
//-------------------------------------- |
//-------------------------------------- |
757 |
void TrkLevel1::Delete(){ |
void TrkLevel1::Delete(){ |
758 |
|
|
759 |
for(Int_t i=0; i<12 ; i++){ |
// Clear(); |
760 |
good[i] = -1; |
if(Cluster)Cluster->Delete(); |
761 |
for(Int_t j=0; j<24 ; j++){ |
if(Cluster)delete Cluster; |
|
cn[j][i] = 0; |
|
|
// cnrms[j][i] = 0; |
|
|
cnn[j][i] = 0; |
|
|
}; |
|
|
}; |
|
|
// |
|
|
Cluster->Delete(); |
|
762 |
|
|
763 |
} |
} |
|
|
|
764 |
//-------------------------------------- |
//-------------------------------------- |
765 |
// |
// |
766 |
// |
// |
767 |
//-------------------------------------- |
//-------------------------------------- |
768 |
TrkCluster *TrkLevel1::GetCluster(int is){ |
TrkCluster *TrkLevel1::GetCluster(int is){ |
769 |
|
|
770 |
// if(is >= this->nclstr()){ |
if(!Cluster)return 0; |
771 |
// cout << "** TrkLevel1::GetCluster(int) ** Cluster "<< is << " does not exits! " << endl; |
if(is >= nclstr()){ |
772 |
// cout << "( Stored clusters nclstr() = "<< this->nclstr()<<" )" << endl; |
cout << "** TrkLevel1::GetCluster(int) ** Cluster "<< is << " does not exits! " << endl; |
773 |
// return 0; |
cout << "( Stored clusters nclstr() = "<< this->nclstr()<<" )" << endl; |
|
// } |
|
|
TClonesArray &t = *(Cluster); |
|
|
TrkCluster *cluster = (TrkCluster*)t[is]; |
|
|
return cluster; |
|
|
} |
|
|
//-------------------------------------- |
|
|
// |
|
|
// |
|
|
//-------------------------------------- |
|
|
/** |
|
|
* Load Position-Finding-Algorythm parameters (call the F77 routine). |
|
|
* |
|
|
*/ |
|
|
int TrkLevel1::LoadPfaParam(TString path){ |
|
|
|
|
|
if( strcmp(path_.path,path.Data()) ){ |
|
|
cout <<"Loading p.f.a. parameters\n"; |
|
|
strcpy(path_.path,path.Data()); |
|
|
path_.pathlen = path.Length(); |
|
|
path_.error = 0; |
|
|
return readetaparam_(); |
|
|
} |
|
774 |
return 0; |
return 0; |
775 |
|
} |
776 |
|
|
777 |
|
TClonesArray &t = *(Cluster); |
778 |
|
TrkCluster *cluster = (TrkCluster*)t[is]; |
779 |
|
return cluster; |
780 |
|
} |
781 |
|
|
782 |
|
|
783 |
|
// int TrkLevel1::GetPfaNbinsAngle(){ |
784 |
|
// TrkParams::Load(4); |
785 |
|
// if( !TrkParams::IsLoaded(4) ){ |
786 |
|
// cout << "int TrkLevel1::GetPfaNbinsAngle() --- ERROR --- p.f.a. parameters not loaded"<<endl; |
787 |
|
// return 0; |
788 |
|
// } |
789 |
|
// return pfa_.nangbin; |
790 |
|
// }; |
791 |
|
|
792 |
|
// int TrkLevel1::GetPfaNbinsETA(){ |
793 |
|
// TrkParams::Load(4); |
794 |
|
// if( !TrkParams::IsLoaded(4) ){ |
795 |
|
// cout << "int TrkLevel1::GetPfaNbinsETA() --- ERROR --- p.f.a. parameters not loaded"<<endl; |
796 |
|
// return 0; |
797 |
|
// } |
798 |
|
// return pfa_.netaval; |
799 |
|
// }; |
800 |
|
|
801 |
|
// /** |
802 |
|
// * |
803 |
|
// * |
804 |
|
// */ |
805 |
|
// float* TrkLevel1::GetPfaCoord(TString pfa, int nview, int nladder, int nang){ |
806 |
|
|
807 |
|
// TrkParams::Load(4); |
808 |
|
// if( !TrkParams::IsLoaded(4) ){ |
809 |
|
// cout << "float* TrkLevel1::GetPfaCoord(TString pfa, int nview, int nladder, int nang) --- ERROR --- p.f.a. parameters not loaded"<<endl; |
810 |
|
// return 0; |
811 |
|
// } |
812 |
|
|
813 |
|
// int nbins = GetPfaNbinsETA(); |
814 |
|
// if(!nbins)return 0; |
815 |
|
|
816 |
|
// float *fcorr = new float [nbins]; |
817 |
|
|
818 |
|
// if(!pfa.CompareTo("ETA2",TString::kIgnoreCase)){ |
819 |
|
// for(int ib=0; ib<nbins; ib++){ |
820 |
|
// fcorr[ib] = pfa_.feta2[nang][nladder][nview][ib]; |
821 |
|
// cout << pfa_.eta2[nang][ib] << " - " << pfa_.feta2[nang][nladder][nview][ib]<<endl;; |
822 |
|
// } |
823 |
|
// }else if (!pfa.CompareTo("ETA3",TString::kIgnoreCase)){ |
824 |
|
// for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.feta3[nang][nladder][nview][ib]; |
825 |
|
// }else if (!pfa.CompareTo("ETA4",TString::kIgnoreCase)){ |
826 |
|
// for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.feta4[nang][nladder][nview][ib]; |
827 |
|
// }else{ |
828 |
|
// cout << pfa<<" pfa parameters not implemented "<<endl; |
829 |
|
// return 0; |
830 |
|
// } |
831 |
|
|
832 |
|
// return fcorr; |
833 |
|
|
834 |
|
// }; |
835 |
|
|
836 |
|
// float* TrkLevel1::GetPfaAbs(TString pfa, int nang){ |
837 |
|
|
838 |
|
// TrkParams::Load(4); |
839 |
|
// if( !TrkParams::IsLoaded(4) ){ |
840 |
|
// cout << "float* TrkLevel1::GetPfaAbs(TString pfa, int nang) --- ERROR --- p.f.a. parameters not loaded"<<endl; |
841 |
|
// return 0; |
842 |
|
// } |
843 |
|
|
844 |
|
// int nbins = GetPfaNbinsETA(); |
845 |
|
// if(!nbins)return 0; |
846 |
|
|
847 |
|
// float *fcorr = new float [nbins]; |
848 |
|
|
849 |
|
// if(!pfa.CompareTo("ETA2",TString::kIgnoreCase)){ |
850 |
|
// for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta2[nang][ib]; |
851 |
|
// }else if (!pfa.CompareTo("ETA3",TString::kIgnoreCase)){ |
852 |
|
// for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta3[nang][ib]; |
853 |
|
// }else if (!pfa.CompareTo("ETA4",TString::kIgnoreCase)){ |
854 |
|
// for(int ib=0; ib<nbins; ib++)fcorr[ib] = pfa_.eta4[nang][ib]; |
855 |
|
// }else{ |
856 |
|
// cout << pfa<<" pfa parameters not implemented "<<endl; |
857 |
|
// return 0; |
858 |
|
// } |
859 |
|
|
860 |
|
// return fcorr; |
861 |
|
|
862 |
|
// }; |
863 |
|
|
864 |
|
/** |
865 |
|
* Method to call the F77 routine that performs level1->level2 processing. |
866 |
|
* The level2 output is stored in a common block, which can be retrieved |
867 |
|
* by mean of the method TrkLevel2::SetFromLevel2Struct(). |
868 |
|
* NB If the TrkLevel1 object is readout from a tree, and the |
869 |
|
* TrkLevel1::ProcessEvent(int pfa) is used to reprocess the event, attention |
870 |
|
* should be payed to the fact that single clusters (clusters not associated |
871 |
|
* with any track) might not be stored. Full reprocessing should be done starting |
872 |
|
* from level0 data. |
873 |
|
*/ |
874 |
|
//int TrkLevel1::ProcessEvent(int pfa){ |
875 |
|
int TrkLevel1::ProcessEvent(){ |
876 |
|
|
877 |
|
// cout << "int TrkLevel1::ProcessEvent()" << endl; |
878 |
|
TrkParams::Load( ); |
879 |
|
if( !TrkParams::IsLoaded() )return 0; |
880 |
|
|
881 |
|
GetLevel1Struct(); |
882 |
|
|
883 |
|
// analysisflight_(&pfa); |
884 |
|
// TrkParams::SetPFA(pfa); |
885 |
|
analysisflight_(); |
886 |
|
|
887 |
|
return 1; |
888 |
|
|
889 |
} |
} |
890 |
|
|
891 |
|
|