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