/[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.10 - (show annotations) (download)
Sat Dec 2 10:42:53 2006 UTC (18 years ago) by pam-fi
Branch: MAIN
Changes since 1.9: +62 -36 lines
implemented a reduced level1 output

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

  ViewVC Help
Powered by ViewVC 1.1.23