Correct memory leak in Find1DClusters
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSPD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16
17 #include "AliITSClusterFinderSPD.h"
18 #include "AliITSMapA1.h"
19 #include "AliITS.h"
20 #include "AliITSdigit.h"
21 #include "AliITSRawCluster.h"
22 #include "AliITSRecPoint.h"
23 #include "AliITSsegmentation.h"
24 #include "AliITSresponse.h"
25 #include "AliRun.h"
26
27
28
29 ClassImp(AliITSClusterFinderSPD)
30
31 //----------------------------------------------------------
32 AliITSClusterFinderSPD::AliITSClusterFinderSPD
33 (AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
34 {
35   // constructor
36     fSegmentation=seg;
37     fDigits=digits;
38     fClusters=recp;
39     fNclusters= fClusters->GetEntriesFast();
40     SetDx();
41     SetDz();
42     fMap=new AliITSMapA1(fSegmentation,fDigits);
43     SetNCells();
44 }
45
46 //_____________________________________________________________________________
47 AliITSClusterFinderSPD::AliITSClusterFinderSPD()
48 {
49   // constructor
50   fSegmentation=0;
51   fDigits=0;
52   fClusters=0;
53   fNclusters=0;
54   fMap=0;
55   SetDx();
56   SetDz();
57   SetNCells();
58   
59 }
60
61 //_____________________________________________________________________________
62 AliITSClusterFinderSPD::~AliITSClusterFinderSPD()
63 {
64   // destructor
65   if (fMap) delete fMap;
66
67
68 }
69 //__________________________________________________________________________
70 AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source){
71   //     Copy Constructor 
72   if(&source == this) return;
73   this->fClusters = source.fClusters ;
74   this->fNclusters = source.fNclusters ;
75   this->fMap = source.fMap ;
76   this->fDz = source.fDz ;
77   this->fDx = source.fDx ;
78   this->fMinNCells = source.fMinNCells ;
79   return;
80 }
81
82 //_________________________________________________________________________
83 AliITSClusterFinderSPD& 
84   AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD &source) {
85   //    Assignment operator
86   if(&source == this) return *this;
87   this->fClusters = source.fClusters ;
88   this->fNclusters = source.fNclusters ;
89   this->fMap = source.fMap ;
90   this->fDz = source.fDz ;
91   this->fDx = source.fDx ;
92   this->fMinNCells = source.fMinNCells ;
93   return *this;
94 }
95
96 //_____________________________________________________________________________
97 void AliITSClusterFinderSPD::SetMap()
98 {
99   // set map
100
101   if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
102
103 }
104
105 //_____________________________________________________________________________
106
107 void AliITSClusterFinderSPD::Find1DClusters()
108 {
109   // Find one dimensional clusters, i.e.
110   // in r*phi(x) direction for each colunm in z direction
111   
112   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
113   
114   // retrieve the parameters 
115   Int_t fNofPixels = fSegmentation->Npz(); 
116   Int_t fMaxNofSamples = fSegmentation->Npx();
117   
118   // read in digits -> do not apply threshold 
119   // signal in fired pixels is always 1
120
121   fMap->FillMap();
122   
123   Int_t nofFoundClusters = 0;
124   
125   Int_t k,it,m;
126   for(k=0;k<fNofPixels;k++) {
127     
128     Int_t mmax = 10;  // a size of the window for the cluster finding
129     
130     for(it=0;it<fMaxNofSamples;it++) {
131       
132       Int_t lclx = 0;
133       Int_t xstart = 0;
134       Int_t xstop = 0;
135       Int_t id = 0;
136       Int_t ilcl =0;
137       
138       for(m=0;m<mmax;m++) {  // find the cluster inside the window
139         id = it+m;
140         if(id >= fMaxNofSamples) break;    // ! no possible for the fadc 
141         
142         if(fMap->TestHit(k,id) == kUnused) {   // start of the cluster
143           lclx += 1;
144           if(lclx == 1) xstart = id;
145           
146         }
147         
148         if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {  
149           // end of cluster if a gap exists
150           xstop = id-1;
151           ilcl = 1;
152           break;
153         }            
154         
155       }   //  end of m-loop
156       
157       if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
158       // continue the "it" loop
159       
160       if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
161       
162       if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {  
163                                    // cluster end is outside of the window,
164         mmax += 5;                 // increase mmax and repeat the cluster
165                                    // finding
166         it -= 1;
167       }
168       
169       if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
170         xstop = fMaxNofSamples - 1;           // the end cluster exists
171         ilcl = 1;
172       } 
173       
174       // ---  Calculate z and x coordinates for one dimensional clusters
175       
176       if(ilcl == 1) {         // new cluster exists
177         it = id;
178         mmax = 10;
179             nofFoundClusters++;
180             Float_t clusterCharge = 0.;
181             Float_t zpitch = fSegmentation->Dpz(k+1); 
182             Float_t clusterZ, dummyX; 
183             Int_t dummy=0;
184             fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ);
185             Float_t zstart = clusterZ - 0.5*zpitch;
186             Float_t zstop = clusterZ + 0.5*zpitch;
187             Float_t clusterX = 0.;
188             Int_t xstartfull = xstart;
189             Int_t xstopfull = xstop;
190             Int_t clusterSizeX = lclx;
191             Int_t clusterSizeZ = 1;
192             
193             Int_t its;
194             for(its=xstart; its<=xstop; its++) {
195               Int_t firedpixel=0;
196               if (fMap->GetHitIndex(k,its)>=0) firedpixel=1; 
197               clusterCharge += firedpixel;
198               clusterX +=its + 0.5;
199             }
200             Float_t fRphiPitch = fSegmentation->Dpx(dummy);
201             clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x 
202             
203             
204             //printf("ClusterZ ClusterX %f %f \n",clusterZ, clusterX);
205             
206             // Int_t nclusters = fClusters->GetEntriesFast();
207             //              cout << nclusters << " clusters" << endl;
208             //            cout<< "Create point"<<endl;
209             
210             // Write the points (coordinates and some cluster information) to the
211             // AliITSRawClusterSPD object
212             
213             AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
214             // fClusters->Add(point);
215             iTS->AddCluster(0,clust);
216             delete clust;
217             //              cout << "Cluster at Ladder: " << fLadder << ", Detector: " <<fDetector<<endl;
218             
219             // cout<<" end of cluster finding for Z pixel "<<endl;
220             
221       }    // new cluster (ilcl=1)
222     } // X direction loop (it)
223   } // Z direction loop (k)
224
225   //fMap->ClearMap();
226   return;
227   
228 }
229
230 //_____________________________________________________________________________
231 void  AliITSClusterFinderSPD::GroupClusters()
232 {
233   // Find two dimensional clusters, i.e. group one dimensional clusters
234   // into two dimensional ones (go both in x and z directions).
235   
236   // get number of clusters for this module
237   Int_t nofClusters = fClusters->GetEntriesFast();
238   nofClusters -= fNclusters;
239   
240   AliITSRawClusterSPD *clusterI;
241   AliITSRawClusterSPD *clusterJ;
242   
243   Int_t *label=new Int_t[nofClusters];  
244   Int_t i,j;
245   for(i=0; i<nofClusters; i++) label[i] = 0;
246   for(i=0; i<nofClusters; i++) {
247     if(label[i] != 0) continue;
248     for(j=i+1; j<nofClusters; j++) { 
249       if(label[j] != 0) continue;
250       clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
251       clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
252       Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
253       if(pair) {     
254         
255         //    if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
256         /*      
257       cout << "clusters " << i << "," << j << " before grouping" << endl;
258               clusterI->PrintInfo();
259               clusterJ->PrintInfo();
260         */    
261         clusterI->Add(clusterJ);
262         //        cout << "remove cluster " << j << endl;
263         label[j] = 1;
264         fClusters->RemoveAt(j);
265         /*
266           cout << "cluster  " << i << " after grouping" << endl;
267           clusterI->PrintInfo();
268         */
269       }  // pair
270     } // J clusters  
271     label[i] = 1;
272   } // I clusters
273   fClusters->Compress();
274   //  Int_t totalNofClusters = fClusters->GetEntriesFast();
275   //  cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
276   
277   delete [] label;
278
279   return;
280   
281   
282 }
283 //_____________________________________________________________________________
284
285 void AliITSClusterFinderSPD::TracksInCluster()
286 {
287   
288   // Find tracks creating one cluster
289
290   // get number of clusters for this module
291   Int_t nofClusters = fClusters->GetEntriesFast();
292   nofClusters -= fNclusters;
293
294   //printf("nofClusters fNclusters %d %d\n",nofClusters,fNclusters);   //commented 18-09-00
295
296   Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz;
297   //  Int_t signal, track0, track1, track2;
298   const Int_t trmax = 100;
299   Int_t cltracks[trmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2;
300
301   for(i=0; i<nofClusters; i++) { 
302     ii = 0;
303     memset(cltracks,-1,sizeof(int)*trmax);
304     tr0=tr1=tr2=-3;
305
306     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
307     //printf("clusterI %p\n",clusterI);               //commented 18-09-00
308
309     nclx = clusterI->NclX();
310     nclz = clusterI->NclZ();
311     xstart = clusterI->XStartf();
312     xstop = clusterI->XStopf();
313     zstart = clusterI->Zend()-nclz+1;
314     zstop = clusterI->Zend();
315
316     Int_t ind; 
317
318      for(iz=0; iz<nclz; iz++) { 
319          jz = zstart + iz;
320        for(ix=0; ix<nclx; ix++) { 
321          jx = xstart + ix;
322          ind = fMap->GetHitIndex(jz,jx);
323          //printf("ind %d\n",ind);                              // commented 18-09-00
324          // this part must be wrong - the index of the digit 
325          // can be zero but not -1 !!! comment out this part
326          // and replace it - see below
327          /*
328          if(ind == 0 && iz >= 0 && ix > 0) {
329           continue;
330          }
331          if(ind == 0 && iz > 0 && ix >= 0) {
332           continue;
333          }
334          if(ind == 0 && iz == 0 && ix == 0 && i > 0) {
335           continue;
336          }
337          */
338          // change the conditions above to :
339          if (ind < 0) continue;
340
341         AliITSdigitSPD *dig = (AliITSdigitSPD*)fMap->GetHit(jz,jx);
342         if (!dig) printf("SPD: something wrong ! dig %p\n",dig);
343         /*
344          signal=dig->fSignal;
345          track0=dig->fTracks[0];
346          track1=dig->fTracks[1];
347          track2=dig->fTracks[2];
348         */
349           for(itr=0; itr<3; itr++) { 
350             tracki = dig->fTracks[itr];
351             //printf("tracki %d\n",tracki);                        //commented 18-09-00
352             if(tracki >= 0) {
353               ii += 1;
354              cltracks[ii-1] = tracki;
355             }
356           }
357        } // ix pixel
358      }  // iz pixel
359  
360      for(is=0; is<trmax; is++) { 
361          if(cltracks[is]<0) continue;
362        for(js=is+1; js<trmax; js++) { 
363          if(cltracks[js]<0) continue;
364          if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
365        }
366      }
367
368      ntr = 0;
369      for(ie=0; ie<trmax; ie++) { 
370        if(cltracks[ie] >= 0) {
371         ntr=ntr+1;
372         if(ntr==1) tr0=cltracks[ie];
373         if(ntr==2) tr1=cltracks[ie];
374         if(ntr==3) tr2=cltracks[ie];
375        }
376      }
377      // if delta ray only
378      if(ntr == 0) ntr = 1;
379
380      clusterI->SetNTracks(ntr);
381      clusterI->SetTracks(tr0,tr1,tr2);
382
383   } // I cluster
384
385 }
386 //_____________________________________________________________________________
387
388 void AliITSClusterFinderSPD::GetRecPoints()
389 {
390   // get rec points
391   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
392   
393   // get number of clusters for this module
394   Int_t nofClusters = fClusters->GetEntriesFast();
395   nofClusters -= fNclusters;
396   const Float_t kconv = 1.0e-4;
397   const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
398   const Float_t kRMSz = 120.0*kconv; // resolution for 425 micron pixel
399
400   Float_t spdLength = fSegmentation->Dz();
401   Float_t spdWidth = fSegmentation->Dx();
402
403   Int_t i;
404   Int_t track0, track1, track2;
405
406   for(i=0; i<nofClusters; i++) { 
407
408     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
409     clusterI->GetTracks(track0, track1, track2); 
410     AliITSRecPoint rnew;
411
412     rnew.SetX((clusterI->X() - spdWidth/2)*kconv);
413     rnew.SetZ((clusterI->Z() - spdLength/2)*kconv);
414     rnew.SetQ(1.);
415     rnew.SetdEdX(0.);
416     rnew.SetSigmaX2(kRMSx*kRMSx);
417     rnew.SetSigmaZ2(kRMSz*kRMSz);
418     rnew.fTracks[0]=track0;
419     rnew.fTracks[1]=track1;
420     rnew.fTracks[2]=track2;
421     iTS->AddRecPoint(rnew);
422   } // I clusters
423
424   fMap->ClearMap();
425   
426 }
427 //_____________________________________________________________________________
428
429 void AliITSClusterFinderSPD::FindRawClusters()
430 {
431   // find raw clusters
432   Find1DClusters();
433   GroupClusters();
434   TracksInCluster();
435   GetRecPoints();
436
437 }
438