New ITS code for new structure and simulations.
[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 "AliRun.h"
19
20
21
22 ClassImp(AliITSClusterFinderSPD)
23
24 //----------------------------------------------------------
25 AliITSClusterFinderSPD::AliITSClusterFinderSPD
26 (AliITSsegmentation *seg, TClonesArray *digits, TClonesArray *recp)   
27 {
28   // constructor
29     fSegmentation=seg;
30     fDigits=digits;
31     fClusters=recp;
32     fNclusters= fClusters->GetEntriesFast();
33     SetDx();
34     SetDz();
35     SetMap();
36     SetNCells();
37 }
38
39 //_____________________________________________________________________________
40 AliITSClusterFinderSPD::AliITSClusterFinderSPD()
41 {
42   // constructor
43   fSegmentation=0;
44   fDigits=0;
45   fClusters=0;
46   fNclusters=0;
47   SetDx();
48   SetDz();
49   SetMap();
50   SetNCells();
51   
52 }
53
54 //__________________________________________________________________________
55 AliITSClusterFinderSPD::AliITSClusterFinderSPD(const AliITSClusterFinderSPD &source){
56   //     Copy Constructor 
57   if(&source == this) return;
58   this->fClusters = source.fClusters ;
59   this->fNclusters = source.fNclusters ;
60   this->fMap = source.fMap ;
61   this->fDz = source.fDz ;
62   this->fDx = source.fDx ;
63   this->fMinNCells = source.fMinNCells ;
64   return;
65 }
66
67 //_________________________________________________________________________
68 AliITSClusterFinderSPD& 
69   AliITSClusterFinderSPD::operator=(const AliITSClusterFinderSPD &source) {
70   //    Assignment operator
71   if(&source == this) return *this;
72   this->fClusters = source.fClusters ;
73   this->fNclusters = source.fNclusters ;
74   this->fMap = source.fMap ;
75   this->fDz = source.fDz ;
76   this->fDx = source.fDx ;
77   this->fMinNCells = source.fMinNCells ;
78   return *this;
79 }
80
81 //_____________________________________________________________________________
82 void AliITSClusterFinderSPD::SetMap()
83 {
84   // set map
85   if(!fMap) fMap=new AliITSMapA1(fSegmentation,fDigits);
86   
87 }
88
89 //_____________________________________________________________________________
90
91 void AliITSClusterFinderSPD::Find1DClusters()
92 {
93   // Find one dimensional clusters, i.e.
94   // in r*phi(x) direction for each colunm in z direction
95   
96   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
97   
98   // retrieve the parameters 
99   Int_t fNofPixels = fSegmentation->Npz(); 
100   Int_t fMaxNofSamples = fSegmentation->Npx();
101   
102   // read in digits -> do not apply threshold 
103   // signal in fired pixels is always 1
104   fMap->FillMap();
105   
106   Int_t nofFoundClusters = 0;
107   
108   Int_t k,it,m;
109   for(k=0;k<fNofPixels;k++) {
110     
111     Int_t mmax = 10;  // a size of the window for the cluster finding
112     
113     for(it=0;it<fMaxNofSamples;it++) {
114       
115       Int_t lclx = 0;
116       Int_t xstart = 0;
117       Int_t xstop = 0;
118       Int_t id = 0;
119       Int_t ilcl =0;
120       
121       for(m=0;m<mmax;m++) {  // find the cluster inside the window
122         id = it+m;
123         if(id >= fMaxNofSamples) break;    // ! no possible for the fadc 
124         
125         if(fMap->GetHitIndex(k,id)) {   // start of the cluster
126           lclx += 1;
127           if(lclx == 1) xstart = id;
128           
129         }
130         
131         if(lclx > 0 && !fMap->GetHitIndex(k,id)) {  
132           // end of cluster if a gap exists
133           xstop = id-1;
134           ilcl = 1;
135           break;
136         }            
137         
138       }   //  end of m-loop
139       
140       if(lclx == 0 && ilcl == 0) it = id; // no cluster in the window,
141       // continue the "it" loop
142       
143       if(id >= fMaxNofSamples && lclx == 0) break; // the x row finished
144       
145       if(id < fMaxNofSamples && ilcl == 0 && lclx > 0) {  
146         // cluster end is outside of the window,
147         mmax += 5;                 // increase mmax and repeat the cluster
148         // finding
149         it -= 1;
150       }
151       
152       if(id >= fMaxNofSamples && lclx > 0) {  // the x row finished but
153         xstop = fMaxNofSamples - 1;          // the end cluster exists
154         ilcl = 1;
155       } 
156       
157       // ---  Calculate z and x coordinates for one dimensional clusters
158       
159       if(ilcl == 1) {         // new cluster exists
160         it = id;
161         mmax = 10;
162             nofFoundClusters++;
163             Float_t clusterCharge = 0.;
164             // get this from segmentation when this will be implemented
165             Float_t zpitch = fSegmentation->Dpz(k+1); 
166             Float_t clusterZ, dummyX; 
167             Int_t dummy=0;
168             fSegmentation->GetCellCxz(dummy,k,dummyX,clusterZ);
169             Float_t zstart = clusterZ - 0.5*zpitch;
170             Float_t zstop = clusterZ + 0.5*zpitch;
171             Float_t clusterX = 0.;
172             Int_t xstartfull = xstart;
173             Int_t xstopfull = xstop;
174             Int_t clusterSizeX = lclx;
175             Int_t clusterSizeZ = 1;
176             
177             Int_t its;
178             for(its=xstart; its<=xstop; its++) {
179               Int_t firedpixel=0;
180               if (fMap->GetHitIndex(k,its)) firedpixel=1; 
181               clusterCharge += firedpixel;
182               clusterX +=its + 0.5;
183             }
184             Float_t fRphiPitch = fSegmentation->Dpx(dummy);
185             clusterX /= (clusterSizeX/fRphiPitch); // center of gravity for x 
186             
187             
188             //printf("ClusterZ ClusterX %f %f \n",clusterZ, clusterX);
189             
190             // Int_t nclusters = fClusters->GetEntriesFast();
191             //              cout << nclusters << " clusters" << endl;
192             //            cout<< "Create point"<<endl;
193             
194             // Write the points (coordinates and some cluster information) to the
195             // AliITSRawClusterSPD object
196             
197             AliITSRawClusterSPD *clust = new AliITSRawClusterSPD(clusterZ,clusterX,clusterCharge,clusterSizeZ,clusterSizeX,xstart,xstop,xstartfull,xstopfull,zstart,zstop,k);
198             // fClusters->Add(point);
199             iTS->AddCluster(0,clust);
200             //              cout << "Cluster at Ladder: " << fLadder << ", Detector: " <<fDetector<<endl;
201             
202             // cout<<" end of cluster finding for Z pixel "<<endl;
203             
204       }    // new cluster (ilcl=1)
205     } // X direction loop (it)
206   } // Z direction loop (k)
207   return;
208   
209 }
210
211 //_____________________________________________________________________________
212 void  AliITSClusterFinderSPD::GroupClusters()
213 {
214   // Find two dimensional clusters, i.e. group one dimensional clusters
215   // into two dimensional ones (go both in x and z directions).
216   
217   // get number of clusters for this module
218   Int_t nofClusters = fClusters->GetEntriesFast();
219   nofClusters -= fNclusters;
220   //printf("Group: fNclusters nofClusters %d %d\n",fNclusters, nofClusters);
221   
222   AliITSRawClusterSPD *clusterI;
223   AliITSRawClusterSPD *clusterJ;
224   
225   //Int_t *label=new Int_t[nofClusters];  // activate this for DEC machines
226   Int_t label[nofClusters];   
227   Int_t i,j;
228   for(i=0; i<nofClusters; i++) label[i] = 0;
229   for(i=0; i<nofClusters; i++) {
230     if(label[i] != 0) continue;
231     for(j=i+1; j<nofClusters; j++) { 
232       if(label[j] != 0) continue;
233       clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
234       clusterJ = (AliITSRawClusterSPD*) fClusters->At(j);
235       Bool_t pair = clusterI->Brother(clusterJ,fDz,fDx);
236       if(pair) {     
237         
238         //    if((clusterI->XStop() == clusterJ->XStart()-1)||(clusterI->XStart()==clusterJ->XStop()+1)) cout<<"!! Diagonal cluster"<<endl;
239         /*    
240               cout << "clusters " << i << "," << j << " before grouping" << endl;
241               clusterI->Print();
242               clusterJ->Print();
243         */    
244         clusterI->Add(clusterJ);
245         //        cout << "remove cluster " << j << endl;
246         label[j] = 1;
247         fClusters->RemoveAt(j);
248         /*
249           cout << "cluster  " << i << " after grouping" << endl;
250           clusterI->Print();
251         */
252       }  // pair
253     } // J clusters  
254     label[i] = 1;
255   } // I clusters
256   fClusters->Compress();
257   //Int_t totalNofClusters = fClusters->GetEntriesFast();
258   //cout << " Nomber of clusters at the group end ="<< totalNofClusters<<endl;
259   
260   return;
261   
262   
263 }
264 //_____________________________________________________________________________
265
266 void AliITSClusterFinderSPD::GetRecPoints()
267 {
268   // get rec points
269   AliITS *iTS=(AliITS*)gAlice->GetModule("ITS");
270   
271   // get number of clusters for this module
272   Int_t nofClusters = fClusters->GetEntriesFast();
273   nofClusters -= fNclusters;
274   //printf("GetRecP: fNclusters nofClusters %d %d\n",fNclusters, nofClusters);
275   
276   const Float_t kconv = 1.0e-4;
277   const Float_t kRMSx = 12.0*kconv; // microns -> cm ITS TDR Table 1.3
278   const Float_t kRMSz = 70.0*kconv; // microns -> cm ITS TDR Table 1.3
279
280   Int_t i;
281   for(i=0; i<nofClusters; i++) { 
282     AliITSRawClusterSPD *clusterI = (AliITSRawClusterSPD*) fClusters->At(i);
283     AliITSRecPoint rnew;
284     rnew.SetX(clusterI->X()*kconv);
285     rnew.SetZ(clusterI->Z()*kconv);
286     rnew.SetQ(1.);
287     rnew.SetdEdX(0.);
288     rnew.SetSigmaX2(kRMSx*kRMSx);
289     rnew.SetSigmaZ2(kRMSz*kRMSz);
290     rnew.SetProbability(1.);
291     iTS->AddRecPoint(rnew);
292   } // I clusters
293   
294 }
295 //_____________________________________________________________________________
296
297 void AliITSClusterFinderSPD::FindRawClusters()
298 {
299   // find raw clusters
300   Find1DClusters();
301   GroupClusters();
302   GetRecPoints();
303 }
304