/[PAMELA software]/DarthVader/TrackerLevel2/src/TrkLevel1.cpp
ViewVC logotype

Contents of /DarthVader/TrackerLevel2/src/TrkLevel1.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (show annotations) (download)
Thu Oct 26 16:22:37 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
Changes since 1.5: +162 -128 lines
fitting methods

1 /**
2 * \file TrkLevel1.cpp
3 * \author Elena Vannuccini
4 */
5 #include <TrkLevel1.h>
6 #include <iostream>
7 using namespace std;
8 //......................................
9 // F77 routines
10 //......................................
11 extern "C" {
12
13 // int readetaparam_();
14 float cog_(int*,int*);
15 float pfaeta_(int*,float*);
16 float pfaeta2_(int*,float*);
17 float pfaeta3_(int*,float*);
18 float pfaeta4_(int*,float*);
19
20 }
21 //--------------------------------------
22 //
23 //
24 //--------------------------------------
25 TrkCluster::TrkCluster(){
26
27 view = 0;
28 maxs = 0;
29 indmax = 0;
30
31 CLlength = 0;
32 clsignal = 0;
33 clsigma = 0;
34 cladc = 0;
35 clbad = 0;
36
37 };
38 //--------------------------------------
39 //
40 //
41 //--------------------------------------
42 TrkCluster::~TrkCluster(){
43
44 delete [] clsignal;
45 delete [] clsigma;
46 delete [] cladc;
47 delete [] clbad;
48 };
49 //--------------------------------------
50 //
51 //
52 //--------------------------------------
53 TrkCluster::TrkCluster(const TrkCluster& t){
54
55 view = t.view;
56 maxs = t.maxs;
57 indmax = t.indmax;
58
59 CLlength = t.CLlength;
60 clsignal = new Float_t[CLlength];
61 clsigma = new Float_t[CLlength];
62 cladc = new Int_t[CLlength];
63 clbad = new Bool_t[CLlength];
64 for(Int_t i=0; i<CLlength;i++){
65 clsignal[i] = t.clsignal[i];
66 clsigma[i] = t.clsigma[i];
67 cladc[i] = t.cladc[i];
68 clbad[i] = t.clbad[i];
69 };
70
71 };
72 //--------------------------------------
73 //
74 //
75 //--------------------------------------
76 /**
77 * Evaluate the cluster signal including a maximum number of adjacent
78 * strips, around maxs, having a significant signal.
79 * @param nstrip Maximum number of strips.
80 * @param cut Inclusion cut ( s > cut*sigma ).
81 * If nstrip<=0 only the inclusion cut is used to determine the cluster size.
82 */
83 Float_t TrkCluster::GetSignal(Int_t nstrip, Float_t cut){
84
85 Float_t s = 0;
86
87 if( nstrip<=0 ){
88 // for(Int_t is = 0; is < CLlength; is++){
89 // Float_t scut = cut*clsigma[is];
90 // if(clsignal[is] > scut) s += clsignal[is];
91 // };
92 for(Int_t is = indmax+1; is < CLlength; is++){
93 Float_t scut = cut*clsigma[is];
94 if(clsignal[is] > scut) s += clsignal[is];
95 else break;
96 };
97 for(Int_t is = indmax; is >=0; is--){
98 Float_t scut = cut*clsigma[is];
99 if(clsignal[is] > scut) s += clsignal[is];
100 else break;
101 };
102 return s;
103 };
104
105
106 Int_t il = indmax;
107 Int_t ir = indmax;
108 Int_t inc = 0;
109
110 if( clsignal[indmax] < cut*clsigma[indmax] ) return 0;
111
112 while ( inc < nstrip ){
113 Float_t sl = -100000;
114 Float_t sr = -100000;
115 if( il >= 0 ) sl = clsignal[il];
116 if( ir < CLlength ) sr = clsignal[ir];
117 if( sl == sr && inc == 0 ){
118 s += clsignal[il]; //cout << inc<<" - "<< clsignal[il]<<" "<<s<<endl;
119 il--;
120 ir++;
121 }else if ( sl >= sr && sl > cut*clsigma[il] && inc !=0 ){
122 s += sl;//cout << inc<<" - "<< clsignal[il]<<" "<<s<<endl;
123 il--;
124 }else if ( sl < sr && sr > cut*clsigma[ir] ){
125 s += sr;//cout << inc<<" - " << clsignal[ir]<<" "<<s<<endl;
126 ir++;
127 }else break;
128
129 inc++;
130 }
131 return s;
132 };
133
134 /**
135 including a ( maximum ) fixed number of adjacent strips (with s>0) around the maxs.
136 * @param nstrip Number of strips.
137 */
138 /**
139 * Evaluate the cluster signal-to-noise, as defined by Turchetta, including a maximum number of adjacent strips, around maxs, having a significant signal.
140 * @param nstrip Maximum number of strips.
141 * @param cut Inclusion cut ( s > cut*sigma ).
142 * If nstrip<=0 only the inclusion cut is used to determine the cluster size.
143 */
144 Float_t TrkCluster::GetSignalToNoise(Int_t nstrip, Float_t cut){
145
146 Float_t sn = 0;
147
148 if( nstrip<=0 ){
149 for(Int_t is = indmax+1; is < CLlength; is++){
150 Float_t scut = cut*clsigma[is];
151 if(clsignal[is] > scut) sn += clsignal[is]/clsigma[is];
152 else break;
153 };
154 for(Int_t is = indmax; is >=0; is--){
155 Float_t scut = cut*clsigma[is];
156 if(clsignal[is] > scut) sn += clsignal[is]/clsigma[is];
157 else break;
158 };
159 return sn;
160 };
161
162
163 Int_t il = indmax;
164 Int_t ir = indmax;
165 Int_t inc = 0;
166
167 if( clsignal[indmax] < cut*clsigma[indmax] ) return 0;
168
169 while ( inc < nstrip ){
170 Float_t sl = -100000;
171 Float_t sr = -100000;
172 if( il >= 0 ) sl = clsignal[il];
173 if( ir < CLlength ) sr = clsignal[ir];
174 if( sl == sr && inc == 0 ){
175 sn += clsignal[il]/clsigma[il];
176 il--;
177 ir++;
178 }else if ( sl >= sr && sl > cut*clsigma[il] && inc !=0 ){
179 sn += sl/clsigma[il];
180 il--;
181 }else if ( sl < sr && sr > cut*clsigma[ir] ){
182 sn += sr/clsigma[ir];
183 ir++;
184 }else break;
185
186 inc++;
187 }
188 return sn;
189 };
190 /**
191 * Evaluate the cluster multiplicity.
192 * @param cut Inclusion cut.
193 */
194 Int_t TrkCluster::GetMultiplicity(Float_t cut){
195 Int_t m = 0;
196 for(Int_t is = 0; is < CLlength; is++){
197 Float_t scut = cut*clsigma[is];
198 if(clsignal[is] > scut) m++;
199 };
200 return m;
201 };
202 /**
203 * True if the cluster contains bad strips.
204 * @param nbad Number of strips around the maximum.
205 */
206 Bool_t TrkCluster::IsBad(Int_t nbad){
207
208 /* Float_t max = 0;
209 Int_t imax = 0;
210 for(Int_t is = 0; is < CLlength; is++){
211 if(clsignal[is] > max){
212 max = clsignal[is];
213 imax = is;
214 };
215 };
216
217 Int_t il,ir;
218 il = imax;
219 ir = imax;*/
220
221 Int_t il,ir;
222 il = indmax;
223 ir = indmax;
224 for(Int_t i=1; i<nbad; i++){
225 if (ir == CLlength && il == 0)break;
226 else if (ir == CLlength && il != 0)il--;
227 else if (ir != CLlength && il == 0)ir++;
228 else{
229 if(clsignal[il-1] > clsignal[ir+1])il--;
230 else ir++;
231 }
232 }
233 Int_t isbad = 0;
234 for(Int_t i=il; i<=ir; i++)isbad += clbad[i];
235
236 return ( isbad != nbad );
237 };
238 //--------------------------------------
239 //
240 //
241 //--------------------------------------
242 void TrkCluster::Dump(){
243
244 cout << "----- Cluster" << endl;
245 cout << "View "<<view << " - Ladder "<<GetLadder()<<endl;
246 cout << "Position of maximun "<< maxs <<endl;
247 cout << "Multiplicity "<< GetMultiplicity() <<endl;
248 cout << "Tot signal "<< GetSignal() << " (ADC channels)"<<endl ;
249 cout << "Signal/Noise "<< GetSignalToNoise();
250 cout <<endl<< "Strip signals ";
251 for(Int_t i =0; i<CLlength; i++)cout << " " <<clsignal[i];
252 cout <<endl<< "Strip sigmas ";
253 for(Int_t i =0; i<CLlength; i++)cout << " " <<clsigma[i];
254 cout <<endl<< "Strip ADC ";
255 for(Int_t i =0; i<CLlength; i++)cout << " " <<cladc[i];
256 cout <<endl<< "Strip BAD ";
257 for(Int_t i =0; i<CLlength; i++){
258 if(i==indmax)cout << " *" <<clbad[i]<<"*";
259 else cout << " " <<clbad[i];
260 }
261 cout << endl;
262
263 }
264 //--------------------------------------
265 //
266 //
267 //--------------------------------------
268 /**
269 * Method to fill a level1 struct with only one cluster (done to use F77 p.f.a. routines on a cluster basis).
270 */
271 cTrkLevel1* TrkCluster::GetLevel1Struct(){
272
273 // cTrkLevel1* l1 = new cTrkLevel1;
274
275 cTrkLevel1* l1 = &level1event_ ;
276
277 l1->nclstr1 = 1;
278 l1->view[0] = view;
279 l1->ladder[0] = GetLadder();
280 l1->maxs[0] = maxs;
281 l1->mult[0] = GetMultiplicity();
282 l1->dedx[0] = GetSignal();
283 l1->indstart[0] = 1;
284 l1->indmax[0] = indmax+1;
285 l1->totCLlength = CLlength;
286 for(Int_t i=0; i<CLlength; i++){
287 l1->clsignal[i] = clsignal[i];
288 l1->clsigma[i] = clsigma[i];
289 l1->cladc[i] = cladc[i];
290 l1->clbad[i] = clbad[i];
291 };
292
293 return l1;
294 };
295 //--------------------------------------
296 //
297 //
298 //--------------------------------------
299 /**
300 * Evaluates the Center-Of-Gravity (COG) of the cluster, in strips, relative to the strip with the maximum signal (TrkCluster::maxs).
301 * @param ncog Number of strips to evaluate COG.
302 * If ncog=0, the COG of the cluster is evaluated according to the cluster multiplicity (defined by the inclusion cut).
303 * If ncog>0, the COG is evaluated using ncog strips, even if they have a negative signal (according to G.Landi)
304 */
305 Float_t TrkCluster::GetCOG(Int_t ncog){
306
307 int ic = 1;
308 GetLevel1Struct();
309 return cog_(&ncog,&ic);
310
311 };
312 /**
313 * Evaluates the Center-Of-Gravity (COG) of the cluster, in strips, relative to the strip with the maximum signal (TrkCluster::maxs),
314 * choosing the number of strips according to the angle, as implemented for the eta-algorythm .
315 * @param angle Projected angle in degree.
316 */
317 Float_t TrkCluster::GetCOG(Float_t angle){
318
319 Int_t neta = 0;
320
321 // Float_t eta = GetETA(0,angle);
322 // for(neta=2; neta<10; neta++) if( eta == GetETA(neta,angle) ) break;
323 // if(eta != GetETA(neta,angle) )cout << "Attenzione!! pasticcio "<<endl;
324
325 if( view%2 ){ //Y
326 neta=2;
327 }else{ //X
328 if( abs(angle) <= 10. ){
329 neta = 2;
330 }else if( abs(angle) > 10. && abs(angle) <= 15. ){
331 neta = 3;
332 }else{
333 neta = 4;
334 };
335 };
336
337 return GetCOG(neta);
338
339 };
340 //--------------------------------------
341 //
342 //
343 //--------------------------------------
344 /**
345 * Evaluates the cluster position, in strips, relative to the strip with the maximum signal (TrkCluster::maxs), by applying the non-linear ETA-algorythm.
346 * @param neta Number of strips to evaluate ETA.
347 * @param angle Projected angle between particle track and detector plane.
348 * Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle.
349 */
350 Float_t TrkCluster::GetETA(Int_t neta, float angle){
351
352 // cout << "GetETA(neta,angle) "<< neta << " "<< angle;
353 // LoadPfaParam();
354
355 float ax = angle;
356 int ic = 1;
357 GetLevel1Struct();
358 if(neta == 0) return pfaeta_(&ic,&ax);
359 else if(neta == 2) return pfaeta2_(&ic,&ax);
360 else if(neta == 3) return pfaeta3_(&ic,&ax);
361 else if(neta == 4) return pfaeta4_(&ic,&ax);
362 else cout << "ETA"<<neta<<" not implemented\n";
363 return 0;
364
365 };
366
367 //--------------------------------------
368 //
369 //
370 //--------------------------------------
371 TrkLevel1::TrkLevel1(){
372
373 Cluster = new TClonesArray("TrkCluster");
374
375 for(Int_t i=0; i<12 ; i++){
376 good[i] = -1;
377 for(Int_t j=0; j<24 ; j++){
378 cn[j][i]=0;
379 // cnrms[j][i]=0;
380 cnn[j][i]=0;
381 };
382 };
383 }
384 //--------------------------------------
385 //
386 //
387 //--------------------------------------
388 void TrkLevel1::Dump(){
389
390 cout<<"DSP status: ";
391 for(Int_t i=0; i<12 ; i++)cout<<good[i]<<" ";
392 cout<<endl;
393
394 TClonesArray &t = *Cluster;
395 for(int i=0; i<this->nclstr(); i++) ((TrkCluster *)t[i])->Dump();
396
397 }
398 //--------------------------------------
399 //
400 //
401 //--------------------------------------
402 /**
403 * Fills a TrkLevel1 object with values from a struct cTrkLevel1 (to get data from F77 common).
404 */
405 void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1){
406
407 // *** CLUSTER ***
408 TrkCluster* t_cl = new TrkCluster();
409 TClonesArray &t = *Cluster;
410 for(int i=0; i<l1->nclstr1; i++){
411
412 t_cl->view = l1->view[i];
413 t_cl->maxs = l1->maxs[i];
414 t_cl->indmax = l1->indmax[i] - l1->indstart[i];
415
416 Int_t from = l1->indstart[i] -1;
417 Int_t to = l1->totCLlength ;
418 if(i != l1->nclstr1-1)to = l1->indstart[i+1] -1 ;
419 t_cl->CLlength = to - from ;
420
421 t_cl->clsignal = new Float_t[t_cl->CLlength];
422 t_cl->clsigma = new Float_t[t_cl->CLlength];
423 t_cl->cladc = new Int_t[t_cl->CLlength];
424 t_cl->clbad = new Bool_t[t_cl->CLlength];
425 Int_t index = 0;
426 for(Int_t is = from; is < to; is++ ){
427 t_cl->clsignal[index] = (Float_t) l1->clsignal[is];
428 t_cl->clsigma[index] = (Float_t) l1->clsigma[is];
429 t_cl->cladc[index] = (Int_t) l1->cladc[is];
430 t_cl->clbad[index] = (Bool_t) l1->clbad[is];
431 index++;
432 };
433
434 new(t[i]) TrkCluster(*t_cl);
435 };
436
437 delete t_cl;
438
439 // ****general variables****
440
441 for(Int_t i=0; i<12 ; i++){
442 good[i] = l1->good[i];
443 for(Int_t j=0; j<24 ; j++){
444 cn[j][i] = l1->cnev[j][i];
445 // cnrms[j][i] = l1->cnrmsev[j][i];
446 cnn[j][i] = l1->cnnev[j][i];
447 };
448 };
449
450 }
451 /**
452 * Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common).
453 */
454
455 cTrkLevel1* TrkLevel1::GetLevel1Struct() {
456
457 cTrkLevel1 *l1=0;
458 //
459 for(Int_t i=0; i<12 ; i++){
460 l1->good[i] = good[i];
461 for(Int_t j=0; j<24 ; j++){
462 l1->cnev[j][i] = cn[j][i];
463 // l1->cnrmsev[j][i] = cnrms[j][i];
464 l1->cnnev[j][i] = cnn[j][i];
465 };
466 };
467
468 // *** CLUSTERS ***
469 l1->nclstr1 = Cluster->GetEntries();
470 for(Int_t i=0;i<l1->nclstr1;i++){
471
472 l1->view[i] = ((TrkCluster *)Cluster->At(i))->view;
473 l1->maxs[i] = ((TrkCluster *)Cluster->At(i))->maxs;
474 // COMPLETARE //
475 // COMPLETARE //
476 // COMPLETARE //
477 // COMPLETARE //
478 // COMPLETARE //
479 // COMPLETARE //
480
481 }
482 // COMPLETARE //
483 // COMPLETARE //
484 // COMPLETARE //
485 // COMPLETARE //
486 // COMPLETARE //
487 // COMPLETARE //
488 return l1;
489 }
490 //--------------------------------------
491 //
492 //
493 //--------------------------------------
494 void TrkLevel1::Clear(){
495
496 for(Int_t i=0; i<12 ; i++){
497 good[i] = -1;
498 for(Int_t j=0; j<24 ; j++){
499 cn[j][i] = 0;
500 // cnrms[j][i] = 0;
501 cnn[j][i] = 0;
502 };
503 };
504 //
505 Cluster->Clear();
506
507 }
508 //--------------------------------------
509 //
510 //
511 //--------------------------------------
512 void TrkLevel1::Delete(){
513
514 for(Int_t i=0; i<12 ; i++){
515 good[i] = -1;
516 for(Int_t j=0; j<24 ; j++){
517 cn[j][i] = 0;
518 // cnrms[j][i] = 0;
519 cnn[j][i] = 0;
520 };
521 };
522 //
523 Cluster->Delete();
524
525 }
526
527 //--------------------------------------
528 //
529 //
530 //--------------------------------------
531 TrkCluster *TrkLevel1::GetCluster(int is){
532
533 if(is >= this->nclstr()){
534 cout << "** TrkLevel1::GetCluster(int) ** Cluster "<< is << " does not exits! " << endl;
535 cout << "( Stored clusters nclstr() = "<< this->nclstr()<<" )" << endl;
536 return 0;
537 }
538 TClonesArray &t = *(Cluster);
539 TrkCluster *cluster = (TrkCluster*)t[is];
540 return cluster;
541 }
542 //--------------------------------------
543 //
544 //
545 //--------------------------------------
546 /**
547 * Load Position-Finding-Algorythm parameters (call the F77 routine).
548 *
549 */
550 int TrkLevel1::LoadPfaParam(TString path){
551
552 if( strcmp(path_.path,path.Data()) ){
553 cout <<"Loading p.f.a. parameters\n";
554 strcpy(path_.path,path.Data());
555 path_.pathlen = path.Length();
556 path_.error = 0;
557 return readetaparam_();
558 }
559 return 0;
560 }
561
562
563 ClassImp(TrkLevel1);
564 ClassImp(TrkCluster);

  ViewVC Help
Powered by ViewVC 1.1.23