Added two missing includes to allow macro compilation (thanks to Laurent for remarkin...
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderV2SPD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, 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 //            Implementation of the ITS clusterer V2 class                //
17 //                                                                        //
18 //          Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch            //
19 //          Unfolding switch from AliITSRecoParam: D. Elia, INFN Bari     //
20 //                                                                        //
21 ////////////////////////////////////////////////////////////////////////////
22
23
24 #include "AliITSCalibrationSPD.h"
25 #include "AliITSClusterFinderV2SPD.h"
26 #include "AliITSRecPoint.h"
27 #include "AliITSgeomTGeo.h"
28 #include "AliITSDetTypeRec.h"
29 #include "AliITSReconstructor.h"
30 #include "AliRawReader.h"
31 #include "AliITSRawStreamSPD.h"
32 #include <TClonesArray.h>
33 #include "AliITSdigitSPD.h"
34 #include "AliITSFOSignalsSPD.h"
35
36 ClassImp(AliITSClusterFinderV2SPD)
37
38 //__________________________________________________________________________
39 AliITSClusterFinderV2SPD::AliITSClusterFinderV2SPD(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),
40 fLastSPD1(AliITSgeomTGeo::GetModuleIndex(2,1,1)-1),
41 fNySPD(256),
42 fNzSPD(160),
43 fYpitchSPD(0.0050),
44 fZ1pitchSPD(0.0425),
45 fZ2pitchSPD(0.0625),
46 fHwSPD(0.64),
47 fHlSPD(3.48){
48
49   //Default constructor
50
51   fYSPD[0]=0.5*fYpitchSPD;
52   for (Int_t m=1; m<fNySPD; m++) fYSPD[m]=fYSPD[m-1]+fYpitchSPD; 
53   fZSPD[0]=fZ1pitchSPD;
54   for (Int_t m=1; m<fNzSPD; m++) {
55     Double_t dz=fZ1pitchSPD;
56     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 || 
57         m==127 || m==128) dz=fZ2pitchSPD; 
58     fZSPD[m]=fZSPD[m-1]+dz;
59   }
60   for (Int_t m=0; m<fNzSPD; m++) {
61     Double_t dz=0.5*fZ1pitchSPD;
62     if (m==31 || m==32 || m==63  || m==64  || m==95 || m==96 || 
63         m==127 || m==128) dz=0.5*fZ2pitchSPD; 
64     fZSPD[m]-=dz;
65   }
66
67 }
68 //__________________________________________________________________________
69 void AliITSClusterFinderV2SPD::FindRawClusters(Int_t mod){
70   //Find clusters V2
71   SetModule(mod);
72   FindClustersSPD(fDigits);
73
74 }
75 //__________________________________________________________________________
76 void AliITSClusterFinderV2SPD::RawdataToClusters(AliRawReader* rawReader, TClonesArray** clusters){
77   //------------------------------------------------------------
78   // This function creates ITS clusters from raw data
79   //------------------------------------------------------------
80   rawReader->Reset();
81   AliITSRawStreamSPD inputSPD(rawReader);
82   FindClustersSPD(&inputSPD, clusters);
83
84 }
85 //__________________________________________________________________________
86 Int_t AliITSClusterFinderV2SPD::ClustersSPD(AliBin* bins, TClonesArray* digits,TClonesArray* clusters,Int_t maxBins,Int_t nzbins,Int_t iModule,Bool_t rawdata){
87   
88   //Cluster finder for SPD (from digits and from rawdata)
89
90   static AliITSRecoParam *repa = NULL;
91   if(!repa){
92     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
93     if(!repa){
94       repa = AliITSRecoParam::GetHighFluxParam();
95       AliWarning("Using default AliITSRecoParam class");
96     }
97   }
98   const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(iModule);
99
100    if (repa->GetSPDRemoveNoisyFlag()) {
101     // Loop on noisy pixels and reset them
102     AliITSCalibrationSPD *cal =  
103       (AliITSCalibrationSPD*) fDetTypeRec->GetCalibrationModel(iModule);
104     for(Int_t ipix = 0; ipix<cal->GetNrBad(); ipix++){
105       Int_t row, col;
106       cal->GetBadPixel(ipix,row,col);
107       //PH      printf(" module %d   row %d  col %d \n",iModule,row,col);
108       Int_t index = (row+1) * nzbins + (col+1);
109       
110       bins[index].SetQ(0);
111       bins[index].SetMask(0xFFFFFFFE);
112     }
113   }
114   
115     if (repa->GetSPDRemoveDeadFlag()) {
116     // Loop on dead pixels and reset them
117     AliITSCalibrationSPD *cal =  
118       (AliITSCalibrationSPD*) fDetTypeRec->GetSPDDeadModel(iModule); 
119     if (cal->IsBad()) return 0; // if all ladder is dead, return to save time
120     for(Int_t ipix = 0; ipix<cal->GetNrBad(); ipix++){
121       Int_t row, col;
122       cal->GetBadPixel(ipix,row,col);
123       Int_t index = (row+1) * nzbins + (col+1);
124       bins[index].SetQ(0);
125       bins[index].SetMask(0xFFFFFFFE);
126     }
127   }
128   
129   Int_t nclu=0;
130   for(Int_t iBin =0; iBin < maxBins;iBin++){
131     if(bins[iBin].IsUsed()) continue;
132     Int_t nBins = 0; 
133     Int_t idxBins[200];
134     FindCluster(iBin,nzbins,bins,nBins,idxBins);
135     if (nBins == 200){
136       Error("ClustersSPD","SPD Too big cluster !\n"); 
137       continue;
138     }
139     Int_t milab[10];
140     for(Int_t ilab=0;ilab<10;ilab++){
141       milab[ilab]=-2;
142     }
143     if(rawdata){
144       milab[3]=fNdet[iModule];
145     }
146     Int_t ymin,ymax,zmin,zmax;
147     if(rawdata){
148       ymin = (idxBins[0] / nzbins) - 1;
149       ymax = ymin;
150       zmin = (idxBins[0] % nzbins) - 1;
151       zmax = zmin;
152     }
153     else{
154       AliITSdigitSPD* dig = (AliITSdigitSPD*)digits->UncheckedAt(idxBins[0]);
155       ymin=dig->GetCoord2();
156       ymax=ymin;
157       zmin=dig->GetCoord1();
158       zmax=zmin;
159     }
160     //PH     if(iModule == 24 || iModule == 25)  printf("\n");
161     for (Int_t idx = 0; idx < nBins; idx++) {
162       Int_t iy;
163       Int_t iz; 
164       if(rawdata){
165         iy  = (idxBins[idx] / nzbins) - 1;
166         iz  = (idxBins[idx] % nzbins) - 1;
167       }
168       else{
169         AliITSdigitSPD* dig = (AliITSdigitSPD*)digits->UncheckedAt(idxBins[idx]);
170         iy = dig->GetCoord2();
171         iz = dig->GetCoord1();
172         //if(iModule == 24 || iModule == 25) printf(" ||  iy %d   iz %d  in Module %d \n",iy,iz,iModule);
173       }
174       if (ymin > iy) ymin = iy;
175       if (ymax < iy) ymax = iy;
176       if (zmin > iz) zmin = iz;
177       if (zmax < iz) zmax = iz;
178
179     }
180     if(!rawdata){
181       for(Int_t l=0;l<nBins;l++){
182         AliITSdigitSPD* dig = (AliITSdigitSPD*)digits->UncheckedAt(idxBins[l]);
183         for(Int_t dlab=0;dlab<10;dlab++){
184           Int_t digitlab = (dig->GetTracks())[dlab];
185           if(digitlab<0) continue;
186           AddLabel(milab,digitlab);
187         }
188         if (milab[9]>0) CheckLabels2(milab);
189       }
190       CheckLabels2(milab);
191     }
192     
193     Int_t idy =0; //max 2 clusters
194     if((iModule <= fLastSPD1) &&idy<3) idy=3;
195     if((iModule > fLastSPD1) &&idy<4) idy=4;
196     Int_t idz=3;
197
198     // Switch the unfolding OFF/ON
199     if(!repa->GetUseUnfoldingInClusterFinderSPD()) {
200       idy=ymax-ymin+1;
201       idz=zmax-zmin+1;
202     }
203     
204     for(Int_t iiz=zmin; iiz<=zmax;iiz+=idz){
205       for(Int_t iiy=ymin;iiy<=ymax;iiy+=idy){
206
207         Int_t ndigits=0;
208         Float_t y=0.,z=0.,q=0.;
209         for(Int_t idx=0;idx<nBins;idx++){
210           Int_t iy;
211           Int_t iz; 
212           if(rawdata){
213             iy  = (idxBins[idx] / nzbins)-1;
214             iz  = (idxBins[idx] % nzbins)-1;
215           }
216           else{
217             AliITSdigitSPD* dig = (AliITSdigitSPD*)digits->UncheckedAt(idxBins[idx]);
218             iy = dig->GetCoord2();
219             iz = dig->GetCoord1();
220           }
221           if(zmax-zmin>=idz || ymax-ymin>=idy){
222             if(TMath::Abs(iy-iiy)>0.75*idy) continue;
223             if(TMath::Abs(iz-iiz)>0.75*idz) continue;
224           }
225           ndigits++;
226           Float_t qBin=0.;
227           if(rawdata) qBin = bins[idxBins[idx]].GetQ();
228           if(!rawdata){
229             AliITSdigitSPD* dig = (AliITSdigitSPD*)digits->UncheckedAt(idxBins[idx]);
230             qBin = (Float_t)dig->GetSignal();
231           }
232           y+= qBin * fYSPD[iy];
233           z+= qBin * fZSPD[iz];
234           q+= qBin;     
235         }// for idx
236         if(ndigits==0) continue;
237          
238         y /= q;
239         z /= q;
240         y -= fHwSPD;
241         z -= fHlSPD;
242
243         Float_t hit[5]; //y,z,sigma(y)^2, sigma(z)^2, charge
244         {
245         Double_t loc[3]={y,0.,z},trk[3]={0.,0.,0.};
246         mT2L->MasterToLocal(loc,trk);
247         hit[0]=trk[1];
248         hit[1]=trk[2];
249         }
250         hit[2] = fYpitchSPD*fYpitchSPD/12.;
251         hit[3] = fZ1pitchSPD*fZ1pitchSPD/12.;
252         hit[4] = 1.;
253
254         if(!rawdata) milab[3]=fNdet[iModule];
255         Int_t info[3] = {ymax-ymin+1,zmax-zmin+1,fNlayer[iModule]};
256         if(!rawdata){
257          AliITSRecPoint cl(milab,hit,info);
258          cl.SetType(nBins);
259          fDetTypeRec->AddRecPoint(cl);
260         }
261         else{
262           Int_t label[4]={milab[0],milab[1],milab[2],milab[3]};
263           AliITSRecPoint cl(label, hit,info);
264           cl.SetType(nBins);
265           new (clusters->AddrAt(nclu)) 
266           AliITSRecPoint(cl);
267         } 
268         nclu++;
269       }// for iiy
270     }// for iiz
271   }//end for iBin
272   return nclu;
273   
274 }
275 //__________________________________________________________________________
276 void AliITSClusterFinderV2SPD::FindClustersSPD(AliITSRawStreamSPD* input, 
277                                         TClonesArray** clusters) 
278 {
279   //------------------------------------------------------------
280   // SPD cluster finder for raw data (this method is called once per event)
281   // Now also fills fast-or fired map
282   //------------------------------------------------------------
283   
284   Int_t nClustersSPD = 0;
285   Int_t kNzBins = fNzSPD + 2;
286   Int_t kNyBins = fNySPD + 2;
287   Int_t kMaxBin = kNzBins * kNyBins;
288   AliBin *binsSPD = new AliBin[kMaxBin];
289   AliBin *binsSPDInit = new AliBin[kMaxBin];  
290   AliBin* bins = NULL;
291
292   // read raw data input stream
293   while (kTRUE) {
294     Bool_t next = input->Next();
295     if (!next || input->IsNewModule()) {
296       Int_t iModule = input->GetPrevModuleID();
297
298       // when all data from a module was read, search for clusters
299       if (bins) { 
300         clusters[iModule] = new TClonesArray("AliITSRecPoint");
301         Int_t nClusters = ClustersSPD(bins,0,clusters[iModule],kMaxBin,kNzBins,iModule,kTRUE);
302         nClustersSPD += nClusters;
303         bins = NULL;
304       }
305
306       if (!next) break;
307       bins = binsSPD;
308       memcpy(binsSPD,binsSPDInit,sizeof(AliBin)*kMaxBin);
309     }
310
311     if (next && bins) {
312       // fill the current digit into the bins array
313       Int_t index = (input->GetCoord2()+1) * kNzBins + (input->GetCoord1()+1);
314       bins[index].SetIndex(index);
315       bins[index].SetQ(1);
316       bins[index].SetMask(1);
317     }
318   }
319
320   delete [] binsSPDInit;
321   delete [] binsSPD;
322   
323   // AliDebug(1,Form("found clusters in ITS SPD: %d", nClustersSPD));
324   Info("FindClustersSPD", "found clusters in ITS SPD: %d", nClustersSPD);
325  
326     // Fill the FastOr fired map
327   for (UInt_t eq=0; eq<20; eq++) {
328     for (UInt_t hs=0; hs<6; hs++) {
329       for (UInt_t chip=0; chip<10; chip++) {
330         if (input->GetFastOrSignal(eq,hs,chip)) {
331           fDetTypeRec->SetFastOrFiredMapOnline(eq,hs,chip);
332         }
333       }
334     }
335   }
336   
337 }
338 //__________________________________________________________________________
339 void AliITSClusterFinderV2SPD::FindClustersSPD(TClonesArray *digits) {
340   //------------------------------------------------------------
341   // SPD cluster finder for digits (this method is called for each module)
342   // Now also fills the fast-or fired map
343   //------------------------------------------------------------
344
345
346   Int_t kNzBins = fNzSPD + 2;
347   const Int_t kMAXBIN=kNzBins*(fNySPD+2);
348
349   Int_t ndigits=digits->GetEntriesFast();
350   AliBin *bins=new AliBin[kMAXBIN];
351
352   Int_t idig;
353   AliITSdigitSPD *digit=0;
354   for (idig=0; idig<ndigits; idig++) {
355     digit=(AliITSdigitSPD*)digits->UncheckedAt(idig);
356     Int_t i=digit->GetCoord2()+1;   //y
357     Int_t j=digit->GetCoord1()+1;
358     Int_t index=i*kNzBins+j;
359
360     bins[index].SetIndex(idig);
361     bins[index].SetQ(1);
362     bins[index].SetMask(1);
363   }
364
365    
366   Int_t nClustersSPD = ClustersSPD(bins,digits,0,kMAXBIN,kNzBins,fModule,kFALSE); 
367   delete [] bins;
368
369   AliDebug(1,Form("found clusters in ITS SPD: %d", nClustersSPD));
370   
371     //  Fill the FastOr fired map
372   AliITSFOSignalsSPD* foSignals = fDetTypeRec->GetFOSignals();
373   if (foSignals) {
374     UInt_t eq = AliITSRawStreamSPD::GetOnlineEqIdFromOffline(fModule);
375     UInt_t hs = AliITSRawStreamSPD::GetOnlineHSFromOffline(fModule);
376     for(UInt_t ch=0; ch<5; ch++) {
377       UInt_t chip = AliITSRawStreamSPD::GetOnlineChipFromOffline(fModule,ch*32);
378       if (foSignals->GetSignal(eq,hs,chip)) {
379         fDetTypeRec->SetFastOrFiredMapOnline(eq,hs,chip);
380       }
381     }
382   }
383   
384 }