/[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.9 - (show annotations) (download)
Fri Nov 24 16:44:15 2006 UTC (18 years ago) by pam-fi
Branch: MAIN
CVS Tags: v2r01
Changes since 1.8: +9 -0 lines
*** empty log message ***

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( fabs(angle) <= 10. ){
329 neta = 2;
330 }else if( fabs(angle) > 10. && fabs(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 cout<<"VA1 mask : "<<endl;
394 for(Int_t i=0; i<12 ; i++){
395 for(Int_t ii=0; ii<24 ; ii++){
396 Int_t mask = cnn[ii][i];
397 if(mask>0)mask=1;
398 cout<<mask<<" ";
399 }
400 cout <<endl;
401 }
402
403 TClonesArray &t = *Cluster;
404 for(int i=0; i<this->nclstr(); i++) ((TrkCluster *)t[i])->Dump();
405
406 }
407 //--------------------------------------
408 //
409 //
410 //--------------------------------------
411 /**
412 * Fills a TrkLevel1 object with values from a struct cTrkLevel1 (to get data from F77 common).
413 */
414 void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1){
415
416 // *** CLUSTER ***
417 TrkCluster* t_cl = new TrkCluster();
418 TClonesArray &t = *Cluster;
419 for(int i=0; i<l1->nclstr1; i++){
420
421 t_cl->view = l1->view[i];
422 t_cl->maxs = l1->maxs[i];
423 t_cl->indmax = l1->indmax[i] - l1->indstart[i];
424
425 Int_t from = l1->indstart[i] -1;
426 Int_t to = l1->totCLlength ;
427 if(i != l1->nclstr1-1)to = l1->indstart[i+1] -1 ;
428 t_cl->CLlength = to - from ;
429
430 t_cl->clsignal = new Float_t[t_cl->CLlength];
431 t_cl->clsigma = new Float_t[t_cl->CLlength];
432 t_cl->cladc = new Int_t[t_cl->CLlength];
433 t_cl->clbad = new Bool_t[t_cl->CLlength];
434 Int_t index = 0;
435 for(Int_t is = from; is < to; is++ ){
436 t_cl->clsignal[index] = (Float_t) l1->clsignal[is];
437 t_cl->clsigma[index] = (Float_t) l1->clsigma[is];
438 t_cl->cladc[index] = (Int_t) l1->cladc[is];
439 t_cl->clbad[index] = (Bool_t) l1->clbad[is];
440 index++;
441 };
442
443 new(t[i]) TrkCluster(*t_cl);
444 };
445
446 delete t_cl;
447
448 // ****general variables****
449
450 for(Int_t i=0; i<12 ; i++){
451 good[i] = l1->good[i];
452 for(Int_t j=0; j<24 ; j++){
453 cn[j][i] = l1->cnev[j][i];
454 // cnrms[j][i] = l1->cnrmsev[j][i];
455 cnn[j][i] = l1->cnnev[j][i];
456 };
457 };
458
459 }
460 /**
461 * Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common).
462 */
463
464 cTrkLevel1* TrkLevel1::GetLevel1Struct() {
465
466 cTrkLevel1 *l1=0;
467 //
468 for(Int_t i=0; i<12 ; i++){
469 l1->good[i] = good[i];
470 for(Int_t j=0; j<24 ; j++){
471 l1->cnev[j][i] = cn[j][i];
472 // l1->cnrmsev[j][i] = cnrms[j][i];
473 l1->cnnev[j][i] = cnn[j][i];
474 };
475 };
476
477 // *** CLUSTERS ***
478 l1->nclstr1 = Cluster->GetEntries();
479 for(Int_t i=0;i<l1->nclstr1;i++){
480
481 l1->view[i] = ((TrkCluster *)Cluster->At(i))->view;
482 l1->maxs[i] = ((TrkCluster *)Cluster->At(i))->maxs;
483 // COMPLETARE //
484 // COMPLETARE //
485 // COMPLETARE //
486 // COMPLETARE //
487 // COMPLETARE //
488 // COMPLETARE //
489
490 }
491 // COMPLETARE //
492 // COMPLETARE //
493 // COMPLETARE //
494 // COMPLETARE //
495 // COMPLETARE //
496 // COMPLETARE //
497 return l1;
498 }
499 //--------------------------------------
500 //
501 //
502 //--------------------------------------
503 void TrkLevel1::Clear(){
504
505 for(Int_t i=0; i<12 ; i++){
506 good[i] = -1;
507 for(Int_t j=0; j<24 ; j++){
508 cn[j][i] = 0;
509 // cnrms[j][i] = 0;
510 cnn[j][i] = 0;
511 };
512 };
513 //
514 Cluster->Clear();
515
516 }
517 //--------------------------------------
518 //
519 //
520 //--------------------------------------
521 void TrkLevel1::Delete(){
522
523 for(Int_t i=0; i<12 ; i++){
524 good[i] = -1;
525 for(Int_t j=0; j<24 ; j++){
526 cn[j][i] = 0;
527 // cnrms[j][i] = 0;
528 cnn[j][i] = 0;
529 };
530 };
531 //
532 Cluster->Delete();
533
534 }
535
536 //--------------------------------------
537 //
538 //
539 //--------------------------------------
540 TrkCluster *TrkLevel1::GetCluster(int is){
541
542 if(is >= this->nclstr()){
543 cout << "** TrkLevel1::GetCluster(int) ** Cluster "<< is << " does not exits! " << endl;
544 cout << "( Stored clusters nclstr() = "<< this->nclstr()<<" )" << endl;
545 return 0;
546 }
547 TClonesArray &t = *(Cluster);
548 TrkCluster *cluster = (TrkCluster*)t[is];
549 return cluster;
550 }
551 //--------------------------------------
552 //
553 //
554 //--------------------------------------
555 /**
556 * Load Position-Finding-Algorythm parameters (call the F77 routine).
557 *
558 */
559 int TrkLevel1::LoadPfaParam(TString path){
560
561 if( strcmp(path_.path,path.Data()) ){
562 cout <<"Loading p.f.a. parameters\n";
563 strcpy(path_.path,path.Data());
564 path_.pathlen = path.Length();
565 path_.error = 0;
566 return readetaparam_();
567 }
568 return 0;
569 }
570
571
572 ClassImp(TrkLevel1);
573 ClassImp(TrkCluster);

  ViewVC Help
Powered by ViewVC 1.1.23