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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by pam-fi, Tue Sep 5 15:15:40 2006 UTC revision 1.4 by pam-fi, Wed Oct 11 08:00:50 2006 UTC
# Line 5  Line 5 
5  #include <TrkLevel1.h>  #include <TrkLevel1.h>
6  #include <iostream>  #include <iostream>
7  using namespace std;  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  //  //
# Line 12  using namespace std; Line 25  using namespace std;
25  TrkCluster::TrkCluster(){  TrkCluster::TrkCluster(){
26                    
27          view     = 0;          view     = 0;
         ladder   = 0;  
28          maxs     = 0;          maxs     = 0;
29          mult     = 0;          indmax   = 0;
30          sgnl     = 0;          
         whichtrk = -1;  
31          CLlength = 0;          CLlength = 0;
32          clsignal = 0;          clsignal = 0;
33          clsigma  = 0;          clsigma  = 0;
# Line 42  TrkCluster::~TrkCluster(){ Line 53  TrkCluster::~TrkCluster(){
53  TrkCluster::TrkCluster(const TrkCluster& t){  TrkCluster::TrkCluster(const TrkCluster& t){
54                    
55          view     = t.view;          view     = t.view;
         ladder   = t.ladder;  
56          maxs     = t.maxs;          maxs     = t.maxs;
57          mult     = t.mult;          indmax   = t.indmax;
         sgnl     = t.sgnl;  
         whichtrk = t.whichtrk;  
58                    
59          CLlength = t.CLlength;            CLlength = t.CLlength;  
60          clsignal = new Float_t[CLlength];          clsignal = new Float_t[CLlength];
# Line 65  TrkCluster::TrkCluster(const TrkCluster& Line 73  TrkCluster::TrkCluster(const TrkCluster&
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[il];
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(){  void TrkCluster::Dump(){
217    
218          cout << "----- Cluster" << endl;          cout << "----- Cluster" << endl;
219          cout << "View "<<view << " - Ladder "<<ladder<<endl;          cout << "View "<<view << " - Ladder "<<GetLadder()<<endl;
220          cout << "(Track ID "<<whichtrk<<")"<<endl;          cout << "Position of maximun "<< maxs <<endl;
221          cout << "Position of maximun "<<maxs<<endl;          cout << "Multiplicity        "<< GetMultiplicity() <<endl;
222          cout << "Multiplicity        "<<mult<<endl;          cout << "Tot signal          "<< GetSignal() << " (ADC channels)"<<endl ;
223          cout << "Tot signal          "<<sgnl<< " (ADC channels)";          cout << "Signal/Noise        "<< GetSignalToNoise();
224          cout <<endl<< "Strip signals       ";          cout <<endl<< "Strip signals       ";
225          for(Int_t i =0; i<CLlength; i++)cout << " " <<clsignal[i];          for(Int_t i =0; i<CLlength; i++)cout << " " <<clsignal[i];
226          cout <<endl<< "Strip sigmas        ";          cout <<endl<< "Strip sigmas        ";
# Line 80  void TrkCluster::Dump(){ Line 228  void TrkCluster::Dump(){
228          cout <<endl<< "Strip ADC           ";          cout <<endl<< "Strip ADC           ";
229          for(Int_t i =0; i<CLlength; i++)cout << " " <<cladc[i];          for(Int_t i =0; i<CLlength; i++)cout << " " <<cladc[i];
230          cout <<endl<< "Strip BAD           ";          cout <<endl<< "Strip BAD           ";
231          for(Int_t i =0; i<CLlength; i++)cout << " " <<clbad[i];          for(Int_t i =0; i<CLlength; i++){
232                    if(i==indmax)cout << "  *" <<clbad[i]<<"*";
233                    else cout << " " <<clbad[i];
234            }
235          cout << endl;          cout << endl;
236                    
237  }  }
# Line 88  void TrkCluster::Dump(){ Line 239  void TrkCluster::Dump(){
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(){  TrkLevel1::TrkLevel1(){
313            
314          good1 = -1;  //      good1 = -1;
315                    
316          Cluster = new TClonesArray("TrkCluster");          Cluster = new TClonesArray("TrkCluster");
317                    
318          for(Int_t i=0; i<12 ; i++){          for(Int_t i=0; i<12 ; i++){
319  //              crc[i] = -1;  //              crc[i] = -1;
320                    good[i] = -1;
321                  for(Int_t j=0; j<24 ; j++){                  for(Int_t j=0; j<24 ; j++){
322                          cnev[j][i]=0;                          cnev[j][i]=0;
323                          cnnev[j][i]=0;                          cnnev[j][i]=0;
324                  };                  };
325                  fshower[i]=0;  //              fshower[i]=0;
326          };          };
327  }  }
328  //--------------------------------------  //--------------------------------------
# Line 109  TrkLevel1::TrkLevel1(){ Line 331  TrkLevel1::TrkLevel1(){
331  //--------------------------------------  //--------------------------------------
332  void TrkLevel1::Dump(){  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;          TClonesArray &t  = *Cluster;
       
339          for(int i=0; i<this->nclstr(); i++)     ((TrkCluster *)t[i])->Dump();          for(int i=0; i<this->nclstr(); i++)     ((TrkCluster *)t[i])->Dump();
340    
341  }  }
# Line 129  void TrkLevel1::SetFromLevel1Struct(cTrk Line 354  void TrkLevel1::SetFromLevel1Struct(cTrk
354          for(int i=0; i<l1->nclstr1; i++){          for(int i=0; i<l1->nclstr1; i++){
355                                    
356                  t_cl->view     = l1->view[i];                  t_cl->view     = l1->view[i];
                 t_cl->ladder   = l1->ladder[i];  
357                  t_cl->maxs     = l1->maxs[i];                  t_cl->maxs     = l1->maxs[i];
358                  t_cl->mult     = l1->mult[i];                  t_cl->indmax   = l1->indmax[i] - l1->indstart[i];
                 t_cl->sgnl     = l1->dedx[i];  
                 t_cl->whichtrk = l1->whichtrack[i]-1;  
359                                    
360                  Int_t from = l1->indstart[i] -1;                  Int_t from = l1->indstart[i] -1;
361                  Int_t to   = l1->totCLlength ;                  Int_t to   = l1->totCLlength ;
# Line 158  void TrkLevel1::SetFromLevel1Struct(cTrk Line 380  void TrkLevel1::SetFromLevel1Struct(cTrk
380                    
381          delete t_cl;          delete t_cl;
382    
383          //  general variables          //  ****general variables****
384    
         good1 = l1->good1;  
385          for(Int_t i=0; i<12 ; i++){          for(Int_t i=0; i<12 ; i++){
386  //              crc[i] = l1->crc[i];                  good[i] = -1;
387                  for(Int_t j=0; j<24 ; j++){                  for(Int_t j=0; j<24 ; j++){
388                          cnev[j][i]     = l1->cnev[j][i];                          cnev[j][i]     = l1->cnev[j][i];
389                          cnnev[j][i] = l1->cnnev[j][i];                          cnnev[j][i] = l1->cnnev[j][i];
390                  };                  };
                 fshower[i] = l1->fshower[i];  
391          };          };
392                    
393  }  }
# Line 175  void TrkLevel1::SetFromLevel1Struct(cTrk Line 395  void TrkLevel1::SetFromLevel1Struct(cTrk
395   * Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common).   * Fills a struct cTrkLevel1 with values from a TrkLevel1 object (to put data into a F77 common).
396   */   */
397    
398  void TrkLevel1::GetLevel1Struct(cTrkLevel1 *l1) const {  cTrkLevel1* TrkLevel1::GetLevel1Struct() {
399              
400          // ********* completare ********* //          cTrkLevel1 *l1=0;
401          // ********* completare ********* //          //
         // ********* completare ********* //  
         // ********* completare ********* //  
         // ********* completare ********* //  
         // ********* completare ********* //  
 //  general variables  
         l1->good1 = good1;  
402          for(Int_t i=0; i<12 ; i++){          for(Int_t i=0; i<12 ; i++){
403  //              l1->crc[i] = crc[i];                  l1->good[i] = good[i];
404                  for(Int_t j=0; j<24 ; j++){                  for(Int_t j=0; j<24 ; j++){
405                          l1->cnev[j][i]     = cnev[j][i];                          l1->cnev[j][i]  = cnev[j][i];
406                          l1->cnnev[j][i] = cnnev[j][i];                          l1->cnnev[j][i] = cnnev[j][i];
407                  };                  };
                 l1->fshower[i] = fshower[i];  
408          };          };
409                    
410  //  *** CLUSTERS ***  //  *** CLUSTERS ***
# Line 199  void TrkLevel1::GetLevel1Struct(cTrkLeve Line 412  void TrkLevel1::GetLevel1Struct(cTrkLeve
412          for(Int_t i=0;i<l1->nclstr1;i++){          for(Int_t i=0;i<l1->nclstr1;i++){
413    
414                  l1->view[i]     = ((TrkCluster *)Cluster->At(i))->view;                  l1->view[i]     = ((TrkCluster *)Cluster->At(i))->view;
                 l1->ladder[i]   = ((TrkCluster *)Cluster->At(i))->ladder;  
415                  l1->maxs[i]     = ((TrkCluster *)Cluster->At(i))->maxs;                  l1->maxs[i]     = ((TrkCluster *)Cluster->At(i))->maxs;
416                  l1->mult[i]     = ((TrkCluster *)Cluster->At(i))->mult;                  // COMPLETARE //
417                  l1->dedx[i]     = ((TrkCluster *)Cluster->At(i))->sgnl;                  // COMPLETARE //
418                    // COMPLETARE //
419                    // COMPLETARE //
420                    // COMPLETARE //
421                    // COMPLETARE //
422                                    
423          }          }
424                    // COMPLETARE //
425          // ********* completare ********* //          // COMPLETARE //
426            // COMPLETARE //
427            // COMPLETARE //
428            // COMPLETARE //
429            // COMPLETARE //
430            return l1;
431  }  }
432  //--------------------------------------  //--------------------------------------
433  //  //
# Line 215  void TrkLevel1::GetLevel1Struct(cTrkLeve Line 435  void TrkLevel1::GetLevel1Struct(cTrkLeve
435  //--------------------------------------  //--------------------------------------
436  void TrkLevel1::Clear(){  void TrkLevel1::Clear(){
437                    
         good1    = -1;  
438          for(Int_t i=0; i<12 ; i++){          for(Int_t i=0; i<12 ; i++){
439  //              crc[i] = -1;                  good[i] = -1;
440                  for(Int_t j=0; j<24 ; j++){                  for(Int_t j=0; j<24 ; j++){
441                          cnev[j][i]     = 0;                          cnev[j][i]     = 0;
442                          cnnev[j][i] = 0;                          cnnev[j][i] = 0;
443                  };                  };
                 fshower[i] = 0;  
444          };          };
445  //      totCLlength = 0;          //
446          Cluster->Clear();          Cluster->Clear();
447    
448  }  }
# Line 232  void TrkLevel1::Clear(){ Line 450  void TrkLevel1::Clear(){
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){  TrkCluster *TrkLevel1::GetCluster(int is){
472    
473          if(is >= this->nclstr()){          if(is >= this->nclstr()){
# Line 243  TrkCluster *TrkLevel1::GetCluster(int is Line 479  TrkCluster *TrkLevel1::GetCluster(int is
479          TrkCluster *cluster = (TrkCluster*)t[is];          TrkCluster *cluster = (TrkCluster*)t[is];
480          return cluster;          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    
 // //--------------------------------------  
 // //  
 // //  
 // //--------------------------------------  
 // TrkTrackRef::TrkTrackRef(){  
 //      for(int ip=0;ip<6;ip++){  
 //              clx[ip]  = 0;  
 //              cly[ip]  = 0;  
 //      };  
 // };  
 // //--------------------------------------  
 // //  
 // //  
 // //--------------------------------------  
 // TrkTrackRef::TrkTrackRef(const TrkTrackRef& t){  
 //      for(int ip=0;ip<6;ip++){  
 //              clx[ip]  = t.clx[ip];  
 //              cly[ip]  = t.cly[ip];  
 //      };  
 // };  
 // //--------------------------------------  
 // //  
 // //  
 // //--------------------------------------  
 // void TrkTrackRef::Clear(){  
 //      for(int ip=0;ip<6;ip++){  
 //              clx[ip]  = 0;  
 //              cly[ip]  = 0;  
 //      };  
 // };  
 // //--------------------------------------  
 // //  
 // //  
 // //--------------------------------------  
 // TrkLevel2Ref::TrkLevel2Ref(){  
 //      Track    = new TClonesArray("TrkTrackRef");  
 //      SingletX = new TClonesArray("TRef");  
 //      SingletY = new TClonesArray("TRef");  
 // };  
 // //--------------------------------------  
 // //  
 // //  
 // //--------------------------------------  
 // void TrkLevel2Ref::SetFromLevel2Struct(cTrkLevel2 *l2){  
 //        
 //      TrkTrackRef*   t_track   = new TrkTrackRef();  
 //      TRef t_singlet = 0;  
 //        
 //      TClonesArray &t = *Track;  
 //      for(int i=0; i<l2->ntrk; i++){  
 //              for(int ip=0;ip<6;ip++){  
 //                      t_track->clx[ip]  = 0;//<<<puntatore al cluster  
 //                      t_track->cly[ip]  = 0;//<<<puntatore al cluster  
 //              };  
 //              new(t[i]) TrkTrackRef(*t_track);  
 //              t_track->Clear();  
 //      };  
 // //  *** SINGLETS ***  
 //      TClonesArray &sx = *SingletX;  
 //      for(int i=0; i<l2->nclsx; i++){  
 //              t_singlet = 0;//<<<puntatore al cluster  
 //              new(sx[i]) TRef(t_singlet);  
 //      }  
 //      TClonesArray &sy = *SingletY;  
 //      for(int i=0; i<l2->nclsy; i++){  
 //              t_singlet = 0;//<<<puntatore al cluster  
 //              new(sy[i]) TRef(t_singlet);  
 //      };  
 //        
 //      delete t_track;  
 // }  
 // //--------------------------------------  
 // //  
 // //  
 // //--------------------------------------  
 // void TrkLevel2Ref::Clear(){  
 //      Track->Clear();  
 //      SingletX->Clear();  
 //      SingletY->Clear();  
 // }  
502    
503  ClassImp(TrkLevel1);  ClassImp(TrkLevel1);
504  ClassImp(TrkCluster);  ClassImp(TrkCluster);
 // ClassImp(TrkTrackRef);  
 // ClassImp(TrkLevel2Ref);  
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23