/[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.11 - (show annotations) (download)
Thu Jan 11 10:20:58 2007 UTC (17 years, 10 months ago) by pam-fi
Branch: MAIN
Changes since 1.10: +186 -163 lines
memory-leak bugs + other bugs fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23