Release version of ITS code
[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   if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
96   
97 }
98
99 //_____________________________________________________________________________
100
101 void AliITSClusterFinderSPD::Find1DClusters()
102 {
103   // Find one dimensional clusters, i.e.
104   // in r*phi(x) direction for each colunm in z direction
105   
106   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
107   
108   // retrieve the parameters 
109   Int_t fNofPixels = fSegmentation->Npz(); 
110   Int_t fMaxNofSamples = fSegmentation->Npx();
111   
112   // read in digits -> do not apply threshold 
113   // signal in fired pixels is always 1
114   fMap->FillMap();
115   
116   Int_t nofFoundClusters = 0;
117   
118   Int_t k,it,m;
119   for(k=0;k<fNofPixels;k++) {
120     
121     Int_t mmax = 10;  // a size of the window for the cluster finding
122     
123     for(it=0;it<fMaxNofSamples;it++) {
124       
125       Int_t lclx = 0;
126       Int_t xstart = 0;
127       Int_t xstop = 0;
128       Int_t id = 0;
129       Int_t ilcl =0;
130       
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         
135         if(fMap->TestHit(k,id) == kUnused) {   // start of the cluster
136           lclx += 1;
137           if(lclx == 1) xstart = id;
138           
139         }
140         
141         if(lclx > 0 && fMap->TestHit(k,id) == kEmpty) {  
142           // end of cluster if a gap exists
143           xstop = id-1;
144           ilcl = 1;
145           break;
146         }            
147         
148       }   //  end of m-loop
149       
150       if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
151       // continue the "it" loop
152       
153       if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
154       
155       if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {  
156                                    // cluster end is outside of the window,
157         mmax += 5;                 // increase mmax and repeat the cluster
158                                    // finding
159         it -= 1;
160       }
161       
162       if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
163         xstop = fMaxNofSamples - 1;           // the end cluster exists
164         ilcl = 1;
165       } 
166       
167       // ---  Calculate z and x coordinates for one dimensional clusters
168       
169       if(ilcl == 1) {         // new cluster exists
170         it = id;
171         mmax = 10;
172             nofFoundClusters++;
173             Float_t clusterCharge = 0.;
174             Float_t zpitch = fSegmentation->Dpz(k+1); 
175             Float_t clusterZ, dummyX; 
176             Int_t dummy=0;
177             fSegmentation->GetPadCxz(dummy,k,dummyX,clusterZ);
178             Float_t zstart = clusterZ - 0.5*zpitch;
179             Float_t zstop = clusterZ + 0.5*zpitch;
180             Float_t clusterX = 0.;
181             Int_t xstartfull = xstart;
182             Int_t xstopfull = xstop;
183             Int_t clusterSizeX = lclx;
184             Int_t clusterSizeZ = 1;
185             
186             Int_t its;
187             for(its=xstart; its<=xstop; its++) {
188               Int_t firedpixel=0;
189               if (fMap->GetHitIndex(k,its)>=0) firedpixel=1; 
190               clusterCharge += firedpixel;
191               clusterX +=its + 0.5;
192             }
193             Float_t fRphiPitch = fSegmentation->Dpx(dummy);
194             clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x 
195             
196             
197             //printf("ClusterZ ClusterX %f %f \n",clusterZ, clusterX);
198             
199             // Int_t nclusters = fClusters->GetEntriesFast();
200             //              cout << nclusters << " clusters" << endl;
201             //            cout<< "Create point"<<endl;
202             
203             // Write the points (coordinates and some cluster information) to the
204             // AliITSRawClusterSPD object
205             
206             AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
207             // fClusters->Add(point);
208             iTS->AddCluster(0,clust);
209             //              cout << "Cluster at Ladder: " << fLadder << ", Detector: " <<fDetector<<endl;
210             
211             // cout<<" end of cluster finding for Z pixel "<<endl;
212             
213       }    // new cluster (ilcl=1)
214     } // X direction loop (it)
215   } // Z direction loop (k)
216
217   //fMap->ClearMap();
218   return;
219   
220 }
221
222 //_____________________________________________________________________________
223 void  AliITSClusterFinderSPD::GroupClusters()
224 {
225   // Find two dimensional clusters, i.e. group one dimensional clusters
226   // into two dimensional ones (go both in x and z directions).
227   
228   // get number of clusters for this module
229   Int_t nofClusters = fClusters->GetEntriesFast();
230   nofClusters -= fNclusters;
231   
232   AliITSRawClusterSPD *clusterI;
233   AliITSRawClusterSPD *clusterJ;
234   
235   Int_t *label=new Int_t[nofClusters];  
236   Int_t i,j;
237   for(i=0; i<nofClusters; i++) label[i] = 0;
238   for(i=0; i<nofClusters; i++) {
239     if(label[i] != 0) continue;
240     for(j=i+1; j<nofClusters; j++) { 
241       if(label[j] != 0) continue;
242       clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
243       clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
244       Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
245       if(pair) {     
246         
247         //    if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
248         /*    
249               cout << "clusters " << i << "," << j << " before grouping" << endl;
250               clusterI->PrintInfo();
251               clusterJ->PrintInfo();
252         */    
253         clusterI->Add(clusterJ);
254         //        cout << "remove cluster " << j << endl;
255         label[j] = 1;
256         fClusters->RemoveAt(j);
257         /*
258           cout << "cluster  " << i << " after grouping" << endl;
259           clusterI->PrintInfo();
260         */
261       }  // pair
262     } // J clusters  
263     label[i] = 1;
264   } // I clusters
265   fClusters->Compress();
266   //Int_t totalNofClusters = fClusters->GetEntriesFast();
267   //cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
268   
269   delete [] label;
270
271   return;
272   
273   
274 }
275 //_____________________________________________________________________________
276
277 void AliITSClusterFinderSPD::GetRecPoints()
278 {
279   // get rec points
280   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
281   
282   // get number of clusters for this module
283   Int_t nofClusters = fClusters->GetEntriesFast();
284   nofClusters -= fNclusters;
285
286   
287   const Float_t kconv = 1.0e-4;
288   const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
289   const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3
290
291   Float_t spdLength = fSegmentation->Dz();
292   Float_t spdWidth = fSegmentation->Dx();
293
294   Int_t i, ix, iz;
295   for(i=0; i<nofClusters; i++) { 
296     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
297     fSegmentation->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
298     AliITSdigitSPD *dig = (AliITSdigitSPD*)fMap->GetHit(iz-1,ix-1);
299     AliITSRecPoint rnew;
300     rnew.SetX((clusterI->X() - spdWidth/2)*kconv);
301     rnew.SetZ((clusterI->Z() - spdLength/2)*kconv);
302     rnew.SetQ(1.);
303     rnew.SetdEdX(0.);
304     rnew.SetSigmaX2(kRMSx*kRMSx);
305     rnew.SetSigmaZ2(kRMSz*kRMSz);
306     rnew.fTracks[0]=dig->fTracks[0];
307     rnew.fTracks[1]=dig->fTracks[1];
308     rnew.fTracks[2]=dig->fTracks[2];
309     iTS->AddRecPoint(rnew);
310   } // I clusters
311
312   fMap->ClearMap();
313   
314 }
315 //_____________________________________________________________________________
316
317 void AliITSClusterFinderSPD::FindRawClusters()
318 {
319   // find raw clusters
320   Find1DClusters();
321   GroupClusters();
322   GetRecPoints();
323
324 }
325