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