Memory leak corrected (M.Bondila)
[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 $Log$
18 Revision 1.15  2001/05/01 22:37:44  nilsen
19 Added a dummy argument to FindRawClusters. Argument used for SSD version.
20
21 Revision 1.14  2001/03/05 14:48:46  nilsen
22 Fixed a reoccuring bug. Array sizes must be declare const.
23
24 */
25
26 #include <iostream.h>
27 #include "AliITSClusterFinderSPD.h"
28 #include "AliITSMapA1.h"
29 #include "AliITS.h"
30 #include "AliITSdigit.h"
31 #include "AliITSRawCluster.h"
32 #include "AliITSRecPoint.h"
33 #include "AliITSsegmentation.h"
34 #include "AliITSresponse.h"
35 #include "AliRun.h"
36
37
38
39 ClassImp(AliITSClusterFinderSPD)
40
41 //----------------------------------------------------------
42 AliITSClusterFinderSPD::AliITSClusterFinderSPD
43 (AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
44 {
45   // constructor
46     fSegmentation=seg;
47     fDigits=digits;
48     fClusters=recp;
49     fNclusters= 0;
50     fMap=new AliITSMapA1(fSegmentation,fDigits);
51     SetDx();
52     SetDz();
53     SetNCells();
54 }
55
56 //_____________________________________________________________________________
57 AliITSClusterFinderSPD::AliITSClusterFinderSPD()
58 {
59   // constructor
60   fSegmentation=0;
61   fDigits=0;
62   fClusters=0;
63   fNclusters=0;
64   fMap=0;
65   SetDx();
66   SetDz();
67   SetNCells();
68   
69 }
70
71 //_____________________________________________________________________________
72 AliITSClusterFinderSPD::~AliITSClusterFinderSPD()
73 {
74   // destructor
75   if (fMap) delete fMap;
76
77
78 }
79 //_____________________________________________________________________________
80
81 void AliITSClusterFinderSPD::Find1DClusters()
82 {
83   // Find one dimensional clusters, i.e.
84   // in r*phi(x) direction for each colunm in z direction
85   
86   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
87   
88   // retrieve the parameters 
89   Int_t fNofPixels = fSegmentation->Npz(); 
90   Int_t fMaxNofSamples = fSegmentation->Npx();
91   
92   // read in digits -> do not apply threshold 
93   // signal in fired pixels is always 1
94
95   fMap->FillMap();
96   
97   Int_t nofFoundClusters = 0;
98   
99   Int_t k,it,m;
100   for(k=0;k<fNofPixels;k++) {
101     
102     Int_t mmax = 10;  // a size of the window for the cluster finding
103     
104     for(it=0;it<fMaxNofSamples;it++) {
105       
106       Int_t lclx = 0;
107       Int_t xstart = 0;
108       Int_t xstop = 0;
109       Int_t id = 0;
110       Int_t ilcl =0;
111       
112       for(m=0;m<mmax;m++) {  // find the cluster inside the window
113         id = it+m;
114         if(id >= fMaxNofSamples) break;    // ! no possible for the fadc 
115         
116         if(fMap->TestHit(k,id) == kUnused) {   // start of the cluster
117           lclx += 1;
118           if(lclx == 1) xstart = id;
119           
120         }
121         
122         if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {  
123           // end of cluster if a gap exists
124           xstop = id-1;
125           ilcl = 1;
126           break;
127         }            
128         
129       }   //  end of m-loop
130       
131       if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
132       // continue the "it" loop
133       
134       if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
135       
136       if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {  
137                                    // cluster end is outside of the window,
138         mmax += 5;                 // increase mmax and repeat the cluster
139                                    // finding
140         it -= 1;
141       }
142       
143       if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
144         xstop = fMaxNofSamples - 1;           // the end cluster exists
145         ilcl = 1;
146       } 
147       
148       // ---  Calculate z and x coordinates for one dimensional clusters
149       
150       if(ilcl == 1) {         // new cluster exists
151         it = id;
152         mmax = 10;
153             nofFoundClusters++;
154             Float_t clusterCharge = 0.;
155             Float_t zpitch = fSegmentation->Dpz(k+1); 
156             Float_t clusterZ, dummyX; 
157             Int_t dummy=0;
158             fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ);
159             Float_t zstart = clusterZ - 0.5*zpitch;
160             Float_t zstop = clusterZ + 0.5*zpitch;
161             Float_t clusterX = 0.;
162             Int_t xstartfull = xstart;
163             Int_t xstopfull = xstop;
164             Int_t clusterSizeX = lclx;
165             Int_t clusterSizeZ = 1;
166             
167             Int_t its;
168             for(its=xstart; its<=xstop; its++) {
169               Int_t firedpixel=0;
170               if (fMap->GetHitIndex(k,its)>=0) firedpixel=1; 
171               clusterCharge += firedpixel;
172               clusterX +=its + 0.5;
173             }
174             Float_t fRphiPitch = fSegmentation->Dpx(dummy);
175             clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x 
176             
177             // Write the points (coordinates and some cluster information) to the
178             // AliITSRawClusterSPD object
179             
180             AliITSRawClusterSPD clust(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
181
182             iTS->AddCluster(0,&clust);
183
184       }    // new cluster (ilcl=1)
185     } // X direction loop (it)
186   } // Z direction loop (k)
187
188   //fMap->ClearMap();
189   return;
190   
191 }
192
193 //_____________________________________________________________________________
194 void  AliITSClusterFinderSPD::GroupClusters()
195 {
196   // Find two dimensional clusters, i.e. group one dimensional clusters
197   // into two dimensional ones (go both in x and z directions).
198   
199   // get number of clusters for this module
200   Int_t nofClusters = fClusters->GetEntriesFast();
201   nofClusters -= fNclusters;
202
203   AliITSRawClusterSPD *clusterI;
204   AliITSRawClusterSPD *clusterJ;
205   
206   Int_t *label=new Int_t[nofClusters];  
207   Int_t i,j;
208   for(i=0; i<nofClusters; i++) label[i] = 0;
209   for(i=0; i<nofClusters; i++) {
210     if(label[i] != 0) continue;
211     for(j=i+1; j<nofClusters; j++) { 
212       if(label[j] != 0) continue;
213       clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
214       clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
215       Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
216       if(pair) {     
217         
218         //    if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
219
220         /*
221               cout << "clusters " << i << "," << j << " before grouping" << endl;
222               clusterI->PrintInfo();
223               clusterJ->PrintInfo();
224         */
225
226         clusterI->Add(clusterJ);
227         //      cout << "remove cluster " << j << endl;
228         label[j] = 1;
229         fClusters->RemoveAt(j);
230
231         /*
232           cout << "cluster  " << i << " after grouping" << endl;
233           clusterI->PrintInfo();
234         */
235
236       }  // pair
237     } // J clusters  
238     label[i] = 1;
239   } // I clusters
240   fClusters->Compress();
241
242   /*
243     Int_t totalNofClusters = fClusters->GetEntriesFast();
244     cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
245   */
246
247   delete [] label;
248
249   return;
250   
251   
252 }
253 //_____________________________________________________________________________
254
255 void AliITSClusterFinderSPD::TracksInCluster()
256 {
257   
258   // Find tracks creating one cluster
259
260   // get number of clusters for this module
261   Int_t nofClusters = fClusters->GetEntriesFast();
262   nofClusters -= fNclusters;
263
264   Int_t i, ix, iz, jx, jz, xstart, xstop, zstart, zstop, nclx, nclz;
265   const Int_t trmax = 100;
266   Int_t cltracks[trmax], itr, tracki, ii, is, js, ie, ntr, tr0, tr1, tr2;
267
268   for(i=0; i<nofClusters; i++) { 
269     ii = 0;
270     memset(cltracks,-1,sizeof(int)*trmax);
271     tr0=tr1=tr2=-1;
272
273     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
274
275     nclx = clusterI->NclX();
276     nclz = clusterI->NclZ();
277     xstart = clusterI->XStartf();
278     xstop = clusterI->XStopf();
279     zstart = clusterI->Zend()-nclz+1;
280     zstop = clusterI->Zend();
281     Int_t ind; 
282
283      for(iz=0; iz<nclz; iz++) { 
284          jz = zstart + iz;
285        for(ix=0; ix<nclx; ix++) { 
286          jx = xstart + ix;
287          ind = fMap->GetHitIndex(jz,jx);
288          if(ind < 0) {
289           continue;
290          }
291          if(ind == 0 && iz >= 0 && ix > 0) {
292           continue;
293          }
294          if(ind == 0 && iz > 0 && ix >= 0) {
295           continue;
296          }
297          if(ind == 0 && iz == 0 && ix == 0 && i > 0) {
298           continue;
299          }
300
301         AliITSdigitSPD *dig = (AliITSdigitSPD*)fMap->GetHit(jz,jx);
302
303         /*
304          signal=dig->fSignal;
305          track0=dig->fTracks[0];
306          track1=dig->fTracks[1];
307          track2=dig->fTracks[2];
308         */
309
310           for(itr=0; itr<3; itr++) { 
311             tracki = dig->fTracks[itr];
312             if(tracki >= 0) {
313               ii += 1;
314               if(ii > 90) { 
315               }
316               if(ii < 99) cltracks[ii-1] = tracki;
317             }
318           }
319        } // ix pixel
320      }  // iz pixel
321  
322      for(is=0; is<trmax; is++) { 
323          if(cltracks[is]<0) continue;
324        for(js=is+1; js<trmax; js++) { 
325          if(cltracks[js]<0) continue;
326          if(cltracks[js]==cltracks[is]) cltracks[js]=-5;
327        }
328      }
329
330      ntr = 0;
331      for(ie=0; ie<trmax; ie++) { 
332        if(cltracks[ie] >= 0) {
333         ntr=ntr+1;
334         if(ntr==1) tr0=cltracks[ie];
335         if(ntr==2) tr1=cltracks[ie];
336         if(ntr==3) tr2=cltracks[ie];
337        }
338      }
339      // if delta ray only
340      if(ntr == 0) ntr = 1;
341
342      clusterI->SetNTracks(ntr);
343      clusterI->SetTracks(tr0,tr1,tr2);
344
345   } // I cluster
346
347 }
348 //_____________________________________________________________________________
349
350 void AliITSClusterFinderSPD::GetRecPoints()
351 {
352   // get rec points
353   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
354   
355   // get number of clusters for this module
356   Int_t nofClusters = fClusters->GetEntriesFast();
357   nofClusters -= fNclusters;
358   const Float_t kconv = 1.0e-4;
359   const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
360   const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3
361
362   Float_t spdLength = fSegmentation->Dz();
363   Float_t spdWidth = fSegmentation->Dx();
364
365   Int_t i;
366   Int_t track0, track1, track2;
367
368   for(i=0; i<nofClusters; i++) { 
369
370     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
371     clusterI->GetTracks(track0, track1, track2); 
372     AliITSRecPoint rnew;
373
374     rnew.SetX((clusterI->X() - spdWidth/2)*kconv);
375     rnew.SetZ((clusterI->Z() - spdLength/2)*kconv);
376     rnew.SetQ(1.);
377     rnew.SetdEdX(0.);
378     rnew.SetSigmaX2(kRMSx*kRMSx);
379     rnew.SetSigmaZ2(kRMSz*kRMSz);
380     rnew.fTracks[0]=track0;
381     rnew.fTracks[1]=track1;
382     rnew.fTracks[2]=track2;
383     iTS->AddRecPoint(rnew);
384   } // I clusters
385
386   fMap->ClearMap();
387   
388 }
389 //_____________________________________________________________________________
390
391 void AliITSClusterFinderSPD::FindRawClusters(Int_t mod)
392 {
393   // find raw clusters
394   Find1DClusters();
395   GroupClusters();
396   TracksInCluster();
397   GetRecPoints();
398
399 }
400
401
402