/[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.5 - (show annotations) (download)
Thu Oct 12 15:41:03 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
Changes since 1.4: +101 -75 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 l1->nclstr1 = 1;
276 l1->view[0] = view;
277 l1->ladder[0] = GetLadder();
278 l1->maxs[0] = maxs;
279 l1->mult[0] = GetMultiplicity();
280 l1->dedx[0] = GetSignal();
281 l1->indstart[0] = 1;
282 l1->indmax[0] = indmax+1;
283 l1->totCLlength = CLlength;
284 for(Int_t i=0; i<CLlength; i++){
285 l1->clsignal[i] = clsignal[i];
286 l1->clsigma[i] = clsigma[i];
287 l1->cladc[i] = cladc[i];
288 l1->clbad[i] = clbad[i];
289 };
290
291 return l1;
292 };
293 //--------------------------------------
294 //
295 //
296 //--------------------------------------
297 /**
298 * Evaluates the Center-Of-Gravity (COG) of the cluster, in strips, relative to the strip with the maximum signal (TrkCluster::maxs).
299 * @param ncog Number of strips to evaluate COG.
300 * If ncog=0, the COG of the cluster is evaluated according to the cluster multiplicity (defined by the inclusion cut).
301 * If ncog>0, the COG is evaluated using ncog strips, even if they have a negative signal (according to G.Landi)
302 */
303 Float_t TrkCluster::GetCOG(Int_t ncog){
304
305 int ic = 1;
306 level1event_ = *GetLevel1Struct();
307 return cog_(&ncog,&ic);
308
309 };
310 //--------------------------------------
311 //
312 //
313 //--------------------------------------
314 /**
315 * Evaluates the cluster position, in strips, relative to the strip with the maximum signal (TrkCluster::maxs), by applying the non-linear ETA-algorythm.
316 * @param neta Number of strips to evaluate ETA.
317 * @param angle Projected angle between particle track and detector plane.
318 * Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle.
319 */
320 Float_t TrkCluster::GetETA(Int_t neta, float angle){
321
322 // LoadPfaParam();
323 int ic = 1;
324 level1event_ = *GetLevel1Struct();
325 if(neta == 0) return pfaeta_(&ic,&angle);
326 else if(neta == 2) return pfaeta2_(&ic,&angle);
327 else if(neta == 3) return pfaeta3_(&ic,&angle);
328 else if(neta == 4) return pfaeta4_(&ic,&angle);
329 else cout << "ETA"<<neta<<" not implemented\n";
330 return 0;
331
332 };
333
334 //--------------------------------------
335 //
336 //
337 //--------------------------------------
338 TrkLevel1::TrkLevel1(){
339
340 // good1 = -1;
341
342 Cluster = new TClonesArray("TrkCluster");
343
344 for(Int_t i=0; i<12 ; i++){
345 // crc[i] = -1;
346 good[i] = -1;
347 for(Int_t j=0; j<24 ; j++){
348 cnev[j][i]=0;
349 cnnev[j][i]=0;
350 };
351 // fshower[i]=0;
352 };
353 }
354 //--------------------------------------
355 //
356 //
357 //--------------------------------------
358 void TrkLevel1::Dump(){
359
360 cout<<"DSP status: ";
361 for(Int_t i=0; i<12 ; i++)cout<<good[i]<<" ";
362 cout<<endl;
363
364 TClonesArray &t = *Cluster;
365 for(int i=0; i<this->nclstr(); i++) ((TrkCluster *)t[i])->Dump();
366
367 }
368 //--------------------------------------
369 //
370 //
371 //--------------------------------------
372 /**
373 * Fills a TrkLevel1 object with values from a struct cTrkLevel1 (to get data from F77 common).
374 */
375 void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1){
376
377 // *** CLUSTER ***
378 TrkCluster* t_cl = new TrkCluster();
379 TClonesArray &t = *Cluster;
380 for(int i=0; i<l1->nclstr1; i++){
381
382 t_cl->view = l1->view[i];
383 t_cl->maxs = l1->maxs[i];
384 t_cl->indmax = l1->indmax[i] - l1->indstart[i];
385
386 Int_t from = l1->indstart[i] -1;
387 Int_t to = l1->totCLlength ;
388 if(i != l1->nclstr1-1)to = l1->indstart[i+1] -1 ;
389 t_cl->CLlength = to - from ;
390
391 t_cl->clsignal = new Float_t[t_cl->CLlength];
392 t_cl->clsigma = new Float_t[t_cl->CLlength];
393 t_cl->cladc = new Int_t[t_cl->CLlength];
394 t_cl->clbad = new Bool_t[t_cl->CLlength];
395 Int_t index = 0;
396 for(Int_t is = from; is < to; is++ ){
397 t_cl->clsignal[index] = (Float_t) l1->clsignal[is];
398 t_cl->clsigma[index] = (Float_t) l1->clsigma[is];
399 t_cl->cladc[index] = (Int_t) l1->cladc[is];
400 t_cl->clbad[index] = (Bool_t) l1->clbad[is];
401 index++;
402 };
403
404 new(t[i]) TrkCluster(*t_cl);
405 };
406
407 delete t_cl;
408
409 // ****general variables****
410
411 for(Int_t i=0; i<12 ; i++){
412 good[i] = -1;
413 for(Int_t j=0; j<24 ; j++){
414 cnev[j][i] = l1->cnev[j][i];
415 cnnev[j][i] = l1->cnnev[j][i];
416 };
417 };
418
419 }
420 /**
421 * Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common).
422 */
423
424 cTrkLevel1* TrkLevel1::GetLevel1Struct() {
425
426 cTrkLevel1 *l1=0;
427 //
428 for(Int_t i=0; i<12 ; i++){
429 l1->good[i] = good[i];
430 for(Int_t j=0; j<24 ; j++){
431 l1->cnev[j][i] = cnev[j][i];
432 l1->cnnev[j][i] = cnnev[j][i];
433 };
434 };
435
436 // *** CLUSTERS ***
437 l1->nclstr1 = Cluster->GetEntries();
438 for(Int_t i=0;i<l1->nclstr1;i++){
439
440 l1->view[i] = ((TrkCluster *)Cluster->At(i))->view;
441 l1->maxs[i] = ((TrkCluster *)Cluster->At(i))->maxs;
442 // COMPLETARE //
443 // COMPLETARE //
444 // COMPLETARE //
445 // COMPLETARE //
446 // COMPLETARE //
447 // COMPLETARE //
448
449 }
450 // COMPLETARE //
451 // COMPLETARE //
452 // COMPLETARE //
453 // COMPLETARE //
454 // COMPLETARE //
455 // COMPLETARE //
456 return l1;
457 }
458 //--------------------------------------
459 //
460 //
461 //--------------------------------------
462 void TrkLevel1::Clear(){
463
464 for(Int_t i=0; i<12 ; i++){
465 good[i] = -1;
466 for(Int_t j=0; j<24 ; j++){
467 cnev[j][i] = 0;
468 cnnev[j][i] = 0;
469 };
470 };
471 //
472 Cluster->Clear();
473
474 }
475 //--------------------------------------
476 //
477 //
478 //--------------------------------------
479 void TrkLevel1::Delete(){
480
481 for(Int_t i=0; i<12 ; i++){
482 good[i] = -1;
483 for(Int_t j=0; j<24 ; j++){
484 cnev[j][i] = 0;
485 cnnev[j][i] = 0;
486 };
487 };
488 //
489 Cluster->Delete();
490
491 }
492
493 //--------------------------------------
494 //
495 //
496 //--------------------------------------
497 TrkCluster *TrkLevel1::GetCluster(int is){
498
499 if(is >= this->nclstr()){
500 cout << "** TrkLevel1::GetCluster(int) ** Cluster "<< is << " does not exits! " << endl;
501 cout << "( Stored clusters nclstr() = "<< this->nclstr()<<" )" << endl;
502 return 0;
503 }
504 TClonesArray &t = *(Cluster);
505 TrkCluster *cluster = (TrkCluster*)t[is];
506 return cluster;
507 }
508 //--------------------------------------
509 //
510 //
511 //--------------------------------------
512 /**
513 * Load Position-Finding-Algorythm parameters (call the F77 routine).
514 *
515 */
516 int TrkLevel1::LoadPfaParam(TString path){
517
518 if( strcmp(path_.path,path.Data()) ){
519 cout <<"Loading p.f.a. parameters\n";
520 strcpy(path_.path,path.Data());
521 path_.pathlen = path.Length();
522 path_.error = 0;
523 return readetaparam_();
524 }
525 return 0;
526 }
527
528
529 ClassImp(TrkLevel1);
530 ClassImp(TrkCluster);

  ViewVC Help
Powered by ViewVC 1.1.23