/[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.4 - (hide annotations) (download)
Wed Oct 11 08:00:50 2006 UTC (18 years, 1 month ago) by pam-fi
Branch: MAIN
Changes since 1.3: +1 -1 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.3 * Evaluate the cluster signal including all adjacent strip with a significant signal ( s > cut*sigma ).
78 pam-fi 1.2 * @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 pam-fi 1.3 * 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 pam-fi 1.2 * @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 pam-fi 1.3 * 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 pam-fi 1.4 sn += clsignal[il]/clsigma[il];
150 pam-fi 1.3 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 pam-fi 1.2 * 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 pam-fi 1.1 void TrkCluster::Dump(){
217    
218     cout << "----- Cluster" << endl;
219 pam-fi 1.2 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 pam-fi 1.1 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 pam-fi 1.2 for(Int_t i =0; i<CLlength; i++){
232     if(i==indmax)cout << " *" <<clbad[i]<<"*";
233     else cout << " " <<clbad[i];
234     }
235 pam-fi 1.1 cout << endl;
236    
237     }
238     //--------------------------------------
239     //
240     //
241     //--------------------------------------
242 pam-fi 1.3 /**
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 pam-fi 1.1 TrkLevel1::TrkLevel1(){
313    
314 pam-fi 1.2 // good1 = -1;
315 pam-fi 1.1
316     Cluster = new TClonesArray("TrkCluster");
317    
318     for(Int_t i=0; i<12 ; i++){
319     // crc[i] = -1;
320 pam-fi 1.2 good[i] = -1;
321 pam-fi 1.1 for(Int_t j=0; j<24 ; j++){
322     cnev[j][i]=0;
323     cnnev[j][i]=0;
324     };
325 pam-fi 1.2 // fshower[i]=0;
326 pam-fi 1.1 };
327     }
328     //--------------------------------------
329     //
330     //
331     //--------------------------------------
332     void TrkLevel1::Dump(){
333    
334 pam-fi 1.2 cout<<"DSP status: ";
335     for(Int_t i=0; i<12 ; i++)cout<<good[i]<<" ";
336     cout<<endl;
337    
338 pam-fi 1.1 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 pam-fi 1.2 t_cl->indmax = l1->indmax[i] - l1->indstart[i];
359 pam-fi 1.1
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 pam-fi 1.2 // ****general variables****
384 pam-fi 1.1
385     for(Int_t i=0; i<12 ; i++){
386 pam-fi 1.2 good[i] = -1;
387 pam-fi 1.1 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 pam-fi 1.3 cTrkLevel1* TrkLevel1::GetLevel1Struct() {
399    
400     cTrkLevel1 *l1=0;
401     //
402 pam-fi 1.1 for(Int_t i=0; i<12 ; i++){
403 pam-fi 1.3 l1->good[i] = good[i];
404 pam-fi 1.1 for(Int_t j=0; j<24 ; j++){
405 pam-fi 1.3 l1->cnev[j][i] = cnev[j][i];
406 pam-fi 1.1 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 pam-fi 1.3 // COMPLETARE //
417     // COMPLETARE //
418     // COMPLETARE //
419     // COMPLETARE //
420     // COMPLETARE //
421     // COMPLETARE //
422 pam-fi 1.1
423     }
424 pam-fi 1.3 // COMPLETARE //
425     // COMPLETARE //
426     // COMPLETARE //
427     // COMPLETARE //
428     // COMPLETARE //
429     // COMPLETARE //
430     return l1;
431     }
432     //--------------------------------------
433     //
434     //
435     //--------------------------------------
436     void TrkLevel1::Clear(){
437 pam-fi 1.1
438 pam-fi 1.3 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 pam-fi 1.1
448     }
449     //--------------------------------------
450     //
451     //
452     //--------------------------------------
453 pam-fi 1.3 void TrkLevel1::Delete(){
454 pam-fi 1.1
455     for(Int_t i=0; i<12 ; i++){
456 pam-fi 1.2 good[i] = -1;
457 pam-fi 1.1 for(Int_t j=0; j<24 ; j++){
458     cnev[j][i] = 0;
459     cnnev[j][i] = 0;
460     };
461     };
462 pam-fi 1.2 //
463 pam-fi 1.3 Cluster->Delete();
464 pam-fi 1.1
465     }
466 pam-fi 1.3
467 pam-fi 1.1 //--------------------------------------
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 pam-fi 1.3 //--------------------------------------
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 pam-fi 1.1
502    
503     ClassImp(TrkLevel1);
504     ClassImp(TrkCluster);

  ViewVC Help
Powered by ViewVC 1.1.23