/[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.3 - (show annotations) (download)
Wed Oct 11 06:53:01 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
Changes since 1.2: +205 -21 lines
some new 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 all adjacent strip with a significant signal ( s > cut*sigma ).
78 * @param cut Inclusion cut.
79 */
80 Float_t TrkCluster::GetSignal(Float_t cut){
81 Float_t s = 0;
82 for(Int_t is = 0; is < CLlength; is++){
83 Float_t scut = cut*clsigma[is];
84 if(clsignal[is] > scut) s += clsignal[is];
85 };
86 return s;
87 };
88 /**
89 * Evaluate the cluster signal including a ( maximum ) fixed number of adjacent strips (with s>0) around the maxs.
90 * @param nstrip Number of strips.
91 */
92 Float_t TrkCluster::GetSignal(Int_t nstrip){
93
94 Float_t s = 0;
95 Int_t il = indmax;
96 Int_t ir = indmax;
97 Int_t inc = 0;
98
99 while ( inc<nstrip ){
100 Float_t sl = 0;
101 Float_t sr = 0;
102 if( il >= 0 ) sl = clsignal[il];
103 if( ir < CLlength ) sr = clsignal[ir];
104 if( sl == sr && inc == 0 ){
105 s += clsignal[il];
106 il--;
107 ir++;
108 }else if ( sl >= sr && sl>0 && inc !=0){
109 s += sl;
110 il--;
111 }else if ( sl < sr && sr>0 ){
112 s += sr;
113 ir++;
114 }else break;
115
116 inc++;
117 }
118 return s;
119 };
120 /**
121 * Evaluate the cluster signal-to-noise, as defined by Turchetta, including all adjacent strip with a significant signal ( s > cut*sigma ).
122 * @param cut Inclusion cut.
123 */
124 Float_t TrkCluster::GetSignalToNoise(Float_t cut){
125 Float_t sn = 0;
126 for(Int_t is = 0; is < CLlength; is++){
127 Float_t scut = cut*clsigma[is];
128 if(clsignal[is] > scut) sn += clsignal[is]/clsigma[is];
129 };
130 return sn;
131 };
132 /**
133 * Evaluate the cluster signal-to-noise, as defined by Turchetta, including a ( maximum ) fixed number of adjacent strips (with s>0) around the maxs.
134 * @param nstrip Number of strips.
135 */
136 Float_t TrkCluster::GetSignalToNoise(Int_t nstrip){
137
138 Float_t sn = 0;
139 Int_t il = indmax;
140 Int_t ir = indmax;
141 Int_t inc = 0;
142
143 while ( inc<nstrip ){
144 Float_t sl = 0;
145 Float_t sr = 0;
146 if( il >= 0 ) sl = clsignal[il];
147 if( ir < CLlength ) sr = clsignal[ir];
148 if( sl == sr && inc == 0 ){
149 sn += clsignal[il]/clsigma[is];
150 il--;
151 ir++;
152 }else if ( sl >= sr && sl>0 && inc !=0){
153 sn += sl/clsigma[il];
154 il--;
155 }else if ( sl < sr && sr>0 ){
156 sn += sr/clsigma[ir];
157 ir++;
158 }else break;
159
160 inc++;
161 }
162 return sn;
163 };
164 /**
165 * Evaluate the cluster multiplicity.
166 * @param cut Inclusion cut.
167 */
168 Int_t TrkCluster::GetMultiplicity(Float_t cut){
169 Int_t m = 0;
170 for(Int_t is = 0; is < CLlength; is++){
171 Float_t scut = cut*clsigma[is];
172 if(clsignal[is] > scut) m++;
173 };
174 return m;
175 };
176 /**
177 * True if the cluster contains bad strips.
178 * @param nbad Number of strips around the maximum.
179 */
180 Bool_t TrkCluster::IsBad(Int_t nbad){
181
182 /* Float_t max = 0;
183 Int_t imax = 0;
184 for(Int_t is = 0; is < CLlength; is++){
185 if(clsignal[is] > max){
186 max = clsignal[is];
187 imax = is;
188 };
189 };
190
191 Int_t il,ir;
192 il = imax;
193 ir = imax;*/
194
195 Int_t il,ir;
196 il = indmax;
197 ir = indmax;
198 for(Int_t i=1; i<nbad; i++){
199 if (ir == CLlength && il == 0)break;
200 else if (ir == CLlength && il != 0)il--;
201 else if (ir != CLlength && il == 0)ir++;
202 else{
203 if(clsignal[il-1] > clsignal[ir+1])il--;
204 else ir++;
205 }
206 }
207 Int_t isbad = 0;
208 for(Int_t i=il; i<=ir; i++)isbad += clbad[i];
209
210 return ( isbad != nbad );
211 };
212 //--------------------------------------
213 //
214 //
215 //--------------------------------------
216 void TrkCluster::Dump(){
217
218 cout << "----- Cluster" << endl;
219 cout << "View "<<view << " - Ladder "<<GetLadder()<<endl;
220 cout << "Position of maximun "<< maxs <<endl;
221 cout << "Multiplicity "<< GetMultiplicity() <<endl;
222 cout << "Tot signal "<< GetSignal() << " (ADC channels)"<<endl ;
223 cout << "Signal/Noise "<< GetSignalToNoise();
224 cout <<endl<< "Strip signals ";
225 for(Int_t i =0; i<CLlength; i++)cout << " " <<clsignal[i];
226 cout <<endl<< "Strip sigmas ";
227 for(Int_t i =0; i<CLlength; i++)cout << " " <<clsigma[i];
228 cout <<endl<< "Strip ADC ";
229 for(Int_t i =0; i<CLlength; i++)cout << " " <<cladc[i];
230 cout <<endl<< "Strip BAD ";
231 for(Int_t i =0; i<CLlength; i++){
232 if(i==indmax)cout << " *" <<clbad[i]<<"*";
233 else cout << " " <<clbad[i];
234 }
235 cout << endl;
236
237 }
238 //--------------------------------------
239 //
240 //
241 //--------------------------------------
242 /**
243 * Method to fill a level1 struct with only one cluster (done to use F77 p.f.a. routines on a cluster basis).
244 */
245 cTrkLevel1* TrkCluster::GetLevel1Struct(){
246
247 cTrkLevel1* l1 = new cTrkLevel1;
248
249 l1->nclstr1 = 1;
250 l1->view[0] = view;
251 l1->ladder[0] = GetLadder();
252 l1->maxs[0] = maxs;
253 l1->mult[0] = GetMultiplicity();
254 l1->dedx[0] = GetSignal();
255 l1->indstart[0] = 1;
256 l1->indmax[0] = indmax+1;
257 l1->totCLlength = CLlength;
258 for(Int_t i=0; i<CLlength; i++){
259 l1->clsignal[i] = clsignal[i];
260 l1->clsigma[i] = clsigma[i];
261 l1->cladc[i] = cladc[i];
262 l1->clbad[i] = clbad[i];
263 };
264
265 return l1;
266 };
267 //--------------------------------------
268 //
269 //
270 //--------------------------------------
271 /**
272 * Evaluates the Center-Of-Gravity (COG) of the cluster, in strips, relative to the strip with the maximum signal (TrkCluster::maxs).
273 * @param ncog Number of strips to evaluate COG.
274 * If ncog=0, the COG of the cluster is evaluated according to the cluster multiplicity (defined by the inclusion cut).
275 * If ncog>0, the COG is evaluated using ncog strips, even if they have a negative signal (according to G.Landi)
276 */
277 Float_t TrkCluster::GetCOG(Int_t ncog){
278
279 int ic = 1;
280 level1event_ = *GetLevel1Struct();
281 return cog_(&ncog,&ic);
282
283 };
284 //--------------------------------------
285 //
286 //
287 //--------------------------------------
288 /**
289 * Evaluates the cluster position, in strips, relative to the strip with the maximum signal (TrkCluster::maxs), by applying the non-linear ETA-algorythm.
290 * @param neta Number of strips to evaluate ETA.
291 * @param angle Projected angle between particle track and detector plane.
292 * Implemented values of neta are 2,3,4. If neta=0, ETA2, ETA3 and ETA4 are applied according to the angle.
293 */
294 Float_t TrkCluster::GetETA(Int_t neta, float angle){
295
296 // LoadPfaParam();
297 int ic = 1;
298 level1event_ = *GetLevel1Struct();
299 if(neta == 0) return pfaeta_(&ic,&angle);
300 else if(neta == 2) return pfaeta2_(&ic,&angle);
301 else if(neta == 3) return pfaeta3_(&ic,&angle);
302 else if(neta == 4) return pfaeta4_(&ic,&angle);
303 else cout << "ETA"<<neta<<" not implemented\n";
304 return 0;
305
306 };
307
308 //--------------------------------------
309 //
310 //
311 //--------------------------------------
312 TrkLevel1::TrkLevel1(){
313
314 // good1 = -1;
315
316 Cluster = new TClonesArray("TrkCluster");
317
318 for(Int_t i=0; i<12 ; i++){
319 // crc[i] = -1;
320 good[i] = -1;
321 for(Int_t j=0; j<24 ; j++){
322 cnev[j][i]=0;
323 cnnev[j][i]=0;
324 };
325 // fshower[i]=0;
326 };
327 }
328 //--------------------------------------
329 //
330 //
331 //--------------------------------------
332 void TrkLevel1::Dump(){
333
334 cout<<"DSP status: ";
335 for(Int_t i=0; i<12 ; i++)cout<<good[i]<<" ";
336 cout<<endl;
337
338 TClonesArray &t = *Cluster;
339 for(int i=0; i<this->nclstr(); i++) ((TrkCluster *)t[i])->Dump();
340
341 }
342 //--------------------------------------
343 //
344 //
345 //--------------------------------------
346 /**
347 * Fills a TrkLevel1 object with values from a struct cTrkLevel1 (to get data from F77 common).
348 */
349 void TrkLevel1::SetFromLevel1Struct(cTrkLevel1 *l1){
350
351 // *** CLUSTER ***
352 TrkCluster* t_cl = new TrkCluster();
353 TClonesArray &t = *Cluster;
354 for(int i=0; i<l1->nclstr1; i++){
355
356 t_cl->view = l1->view[i];
357 t_cl->maxs = l1->maxs[i];
358 t_cl->indmax = l1->indmax[i] - l1->indstart[i];
359
360 Int_t from = l1->indstart[i] -1;
361 Int_t to = l1->totCLlength ;
362 if(i != l1->nclstr1-1)to = l1->indstart[i+1] -1 ;
363 t_cl->CLlength = to - from ;
364
365 t_cl->clsignal = new Float_t[t_cl->CLlength];
366 t_cl->clsigma = new Float_t[t_cl->CLlength];
367 t_cl->cladc = new Int_t[t_cl->CLlength];
368 t_cl->clbad = new Bool_t[t_cl->CLlength];
369 Int_t index = 0;
370 for(Int_t is = from; is < to; is++ ){
371 t_cl->clsignal[index] = (Float_t) l1->clsignal[is];
372 t_cl->clsigma[index] = (Float_t) l1->clsigma[is];
373 t_cl->cladc[index] = (Int_t) l1->cladc[is];
374 t_cl->clbad[index] = (Bool_t) l1->clbad[is];
375 index++;
376 };
377
378 new(t[i]) TrkCluster(*t_cl);
379 };
380
381 delete t_cl;
382
383 // ****general variables****
384
385 for(Int_t i=0; i<12 ; i++){
386 good[i] = -1;
387 for(Int_t j=0; j<24 ; j++){
388 cnev[j][i] = l1->cnev[j][i];
389 cnnev[j][i] = l1->cnnev[j][i];
390 };
391 };
392
393 }
394 /**
395 * Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common).
396 */
397
398 cTrkLevel1* TrkLevel1::GetLevel1Struct() {
399
400 cTrkLevel1 *l1=0;
401 //
402 for(Int_t i=0; i<12 ; i++){
403 l1->good[i] = good[i];
404 for(Int_t j=0; j<24 ; j++){
405 l1->cnev[j][i] = cnev[j][i];
406 l1->cnnev[j][i] = cnnev[j][i];
407 };
408 };
409
410 // *** CLUSTERS ***
411 l1->nclstr1 = Cluster->GetEntries();
412 for(Int_t i=0;i<l1->nclstr1;i++){
413
414 l1->view[i] = ((TrkCluster *)Cluster->At(i))->view;
415 l1->maxs[i] = ((TrkCluster *)Cluster->At(i))->maxs;
416 // COMPLETARE //
417 // COMPLETARE //
418 // COMPLETARE //
419 // COMPLETARE //
420 // COMPLETARE //
421 // COMPLETARE //
422
423 }
424 // COMPLETARE //
425 // COMPLETARE //
426 // COMPLETARE //
427 // COMPLETARE //
428 // COMPLETARE //
429 // COMPLETARE //
430 return l1;
431 }
432 //--------------------------------------
433 //
434 //
435 //--------------------------------------
436 void TrkLevel1::Clear(){
437
438 for(Int_t i=0; i<12 ; i++){
439 good[i] = -1;
440 for(Int_t j=0; j<24 ; j++){
441 cnev[j][i] = 0;
442 cnnev[j][i] = 0;
443 };
444 };
445 //
446 Cluster->Clear();
447
448 }
449 //--------------------------------------
450 //
451 //
452 //--------------------------------------
453 void TrkLevel1::Delete(){
454
455 for(Int_t i=0; i<12 ; i++){
456 good[i] = -1;
457 for(Int_t j=0; j<24 ; j++){
458 cnev[j][i] = 0;
459 cnnev[j][i] = 0;
460 };
461 };
462 //
463 Cluster->Delete();
464
465 }
466
467 //--------------------------------------
468 //
469 //
470 //--------------------------------------
471 TrkCluster *TrkLevel1::GetCluster(int is){
472
473 if(is >= this->nclstr()){
474 cout << "** TrkLevel1::GetCluster(int) ** Cluster "<< is << " does not exits! " << endl;
475 cout << "( Stored clusters nclstr() = "<< this->nclstr()<<" )" << endl;
476 return 0;
477 }
478 TClonesArray &t = *(Cluster);
479 TrkCluster *cluster = (TrkCluster*)t[is];
480 return cluster;
481 }
482 //--------------------------------------
483 //
484 //
485 //--------------------------------------
486 /**
487 * Load Position-Finding-Algorythm parameters (call the F77 routine).
488 *
489 */
490 int TrkLevel1::LoadPfaParam(TString path){
491
492 if( strcmp(path_.path,path.Data()) ){
493 cout <<"Loading p.f.a. parameters\n";
494 strcpy(path_.path,path.Data());
495 path_.pathlen = path.Length();
496 path_.error = 0;
497 return readetaparam_();
498 }
499 return 0;
500 }
501
502
503 ClassImp(TrkLevel1);
504 ClassImp(TrkCluster);

  ViewVC Help
Powered by ViewVC 1.1.23