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