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