/[PAMELA software]/DarthVader/TrackerLevel2/src/TrkLevel1.cpp
ViewVC logotype

Annotation of /DarthVader/TrackerLevel2/src/TrkLevel1.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (hide annotations) (download)
Thu Oct 12 15:41:03 2006 UTC (18 years, 2 months ago) by pam-fi
Branch: MAIN
Changes since 1.4: +101 -75 lines
*** empty log message ***

1 pam-fi 1.1 /**
2     * \file TrkLevel1.cpp
3     * \author Elena Vannuccini
4     */
5     #include <TrkLevel1.h>
6     #include <iostream>
7     using namespace std;
8 pam-fi 1.3 //......................................
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 pam-fi 1.1 //--------------------------------------
22     //
23     //
24     //--------------------------------------
25     TrkCluster::TrkCluster(){
26    
27     view = 0;
28     maxs = 0;
29 pam-fi 1.2 indmax = 0;
30    
31 pam-fi 1.1 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 pam-fi 1.2 indmax = t.indmax;
58 pam-fi 1.1
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 pam-fi 1.2 /**
77 pam-fi 1.5 * 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 pam-fi 1.2 */
83 pam-fi 1.5 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 pam-fi 1.2 };
102     return s;
103 pam-fi 1.5 };
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 pam-fi 1.2 };
133 pam-fi 1.5
134 pam-fi 1.2 /**
135 pam-fi 1.5 including a ( maximum ) fixed number of adjacent strips (with s>0) around the maxs.
136 pam-fi 1.3 * @param nstrip Number of strips.
137     */
138 pam-fi 1.2 /**
139 pam-fi 1.5 * 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 pam-fi 1.3 */
144 pam-fi 1.5 Float_t TrkCluster::GetSignalToNoise(Int_t nstrip, Float_t cut){
145 pam-fi 1.3
146 pam-fi 1.5 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 pam-fi 1.3 return sn;
160 pam-fi 1.5 };
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 pam-fi 1.3 };
190     /**
191 pam-fi 1.2 * 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 pam-fi 1.1 void TrkCluster::Dump(){
243    
244     cout << "----- Cluster" << endl;
245 pam-fi 1.2 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 pam-fi 1.1 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 pam-fi 1.2 for(Int_t i =0; i<CLlength; i++){
258     if(i==indmax)cout << " *" <<clbad[i]<<"*";
259     else cout << " " <<clbad[i];
260     }
261 pam-fi 1.1 cout << endl;
262    
263     }
264     //--------------------------------------
265     //
266     //
267     //--------------------------------------
268 pam-fi 1.3 /**
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 pam-fi 1.1 TrkLevel1::TrkLevel1(){
339    
340 pam-fi 1.2 // good1 = -1;
341 pam-fi 1.1
342     Cluster = new TClonesArray("TrkCluster");
343    
344     for(Int_t i=0; i<12 ; i++){
345     // crc[i] = -1;
346 pam-fi 1.2 good[i] = -1;
347 pam-fi 1.1 for(Int_t j=0; j<24 ; j++){
348     cnev[j][i]=0;
349     cnnev[j][i]=0;
350     };
351 pam-fi 1.2 // fshower[i]=0;
352 pam-fi 1.1 };
353     }
354     //--------------------------------------
355     //
356     //
357     //--------------------------------------
358     void TrkLevel1::Dump(){
359    
360 pam-fi 1.2 cout<<"DSP status: ";
361     for(Int_t i=0; i<12 ; i++)cout<<good[i]<<" ";
362     cout<<endl;
363    
364 pam-fi 1.1 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 pam-fi 1.2 t_cl->indmax = l1->indmax[i] - l1->indstart[i];
385 pam-fi 1.1
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 pam-fi 1.2 // ****general variables****
410 pam-fi 1.1
411     for(Int_t i=0; i<12 ; i++){
412 pam-fi 1.2 good[i] = -1;
413 pam-fi 1.1 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 pam-fi 1.3 cTrkLevel1* TrkLevel1::GetLevel1Struct() {
425    
426     cTrkLevel1 *l1=0;
427     //
428 pam-fi 1.1 for(Int_t i=0; i<12 ; i++){
429 pam-fi 1.3 l1->good[i] = good[i];
430 pam-fi 1.1 for(Int_t j=0; j<24 ; j++){
431 pam-fi 1.3 l1->cnev[j][i] = cnev[j][i];
432 pam-fi 1.1 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 pam-fi 1.3 // COMPLETARE //
443     // COMPLETARE //
444     // COMPLETARE //
445     // COMPLETARE //
446     // COMPLETARE //
447     // COMPLETARE //
448 pam-fi 1.1
449     }
450 pam-fi 1.3 // COMPLETARE //
451     // COMPLETARE //
452     // COMPLETARE //
453     // COMPLETARE //
454     // COMPLETARE //
455     // COMPLETARE //
456     return l1;
457     }
458     //--------------------------------------
459     //
460     //
461     //--------------------------------------
462     void TrkLevel1::Clear(){
463 pam-fi 1.1
464 pam-fi 1.3 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 pam-fi 1.1
474     }
475     //--------------------------------------
476     //
477     //
478     //--------------------------------------
479 pam-fi 1.3 void TrkLevel1::Delete(){
480 pam-fi 1.1
481     for(Int_t i=0; i<12 ; i++){
482 pam-fi 1.2 good[i] = -1;
483 pam-fi 1.1 for(Int_t j=0; j<24 ; j++){
484     cnev[j][i] = 0;
485     cnnev[j][i] = 0;
486     };
487     };
488 pam-fi 1.2 //
489 pam-fi 1.3 Cluster->Delete();
490 pam-fi 1.1
491     }
492 pam-fi 1.3
493 pam-fi 1.1 //--------------------------------------
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 pam-fi 1.3 //--------------------------------------
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 pam-fi 1.1
528    
529     ClassImp(TrkLevel1);
530     ClassImp(TrkCluster);

  ViewVC Help
Powered by ViewVC 1.1.23