29e9714ebcc356cdff9498445ae8780a0a875855
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSDDfast.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
17
18 ////////////////////////////////////////////////////////////////////////////
19 //         Implementation of the ITS SDD fast clusterer  class            //
20 //                                                                        //
21 //   Origin: Simone Capodicasa, Universita e INFN, capodica@to.infn.it    //
22 //                                                                        //
23 ////////////////////////////////////////////////////////////////////////////
24
25 #include <vector>
26 #include <TClonesArray.h>
27 #include <TBits.h>
28 #include <TMath.h>
29 #include <TH2F.h>
30 #include <TFile.h>
31 #include "AliITSClusterFinderSDDfast.h"
32 #include "AliITSRecPoint.h"
33 #include "AliITSRecPointContainer.h"
34 #include "AliITSDetTypeRec.h"
35 #include "AliRawReader.h"
36 #include "AliITSRawStreamSDD.h"
37 #include "AliITSRawStreamSDDCompressed.h"
38 #include "AliITSCalibrationSDD.h"
39 #include "AliITSresponseSDD.h"
40 #include "AliITSDetTypeRec.h"
41 #include "AliITSReconstructor.h"
42 #include "AliITSsegmentationSDD.h"
43 #include "AliITSdigitSDD.h"
44 #include "AliITSgeomTGeo.h"
45
46 ClassImp(AliITSClusterFinderSDDfast)
47
48 AliITSClusterFinderSDDfast::AliITSClusterFinderSDDfast(AliITSDetTypeRec* dettyp):AliITSClusterFinder(dettyp),
49   fNAnodes(0),
50   fNTimeBins(0),
51   fNZbins(0),
52   fNXbins(0),
53   fDDLBins(),
54   fCutOnPeakLoose(0.),
55   fCutOnPeakTight(0.),
56   fMaxDrTimeForTightCut(0.)
57 {
58   //Default constructor
59   fNAnodes = GetSeg()->NpzHalf();
60   fNZbins = fNAnodes+2;
61   fNTimeBins = GetSeg()->Npx();
62   fNXbins = fNTimeBins+2;
63   AliDebug(2,Form("Cells in SDD cluster finder: Andoes=%d  TimeBins=%d",fNAnodes,fNTimeBins));
64   SetPeakSelection(15.,30.,2000.);
65   fDDLBins.resize(kHybridsPerDDL);
66 }
67
68
69 //______________________________________________________________________
70 AliITSClusterFinderSDDfast::~AliITSClusterFinderSDDfast()
71 {
72   //Default destructor
73 }
74
75 //______________________________________________________________________
76 void AliITSClusterFinderSDDfast::FindRawClusters(Int_t mod){
77
78   //Find clusters 
79   SetModule(mod);
80   FindClustersSDD(fDigits);
81
82 }
83
84 //______________________________________________________________________
85 void AliITSClusterFinderSDDfast::FindClustersSDD(TClonesArray *digits){
86
87   std::vector<int> bins0;
88   std::vector<int> bins1;
89   const Int_t kMapDim=fNZbins*fNXbins/32;
90   Int_t map0[kMapDim];
91   Int_t map1[kMapDim];
92   for(Int_t j=0;j<kMapDim;++j){
93     map0[j]=map1[j]=0;
94   }
95   AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(fModule);
96   if(cal==0){
97     AliError(Form("Calibration object not present for SDD module %d\n",fModule));
98     return;
99   }
100
101   AliITSdigitSDD *d=0;
102   Int_t i, ndigits=digits->GetEntriesFast();
103   for (i=0; i<ndigits; i++){
104     d=(AliITSdigitSDD*)digits->UncheckedAt(i);
105     Int_t ian=d->GetCoord1();
106     Int_t itb=d->GetCoord2();
107     Float_t gain=cal->GetChannelGain(ian)/fDetTypeRec->GetAverageGainSDD();
108     Float_t charge=d->GetSignal(); // returns expanded signal
109     // (10 bit, low threshold already added)
110     Float_t baseline = cal->GetBaseline(ian);
111     if(charge>baseline) charge-=baseline;
112     else charge=0;
113
114     if(gain>0.){ // Bad channels have gain=0.
115       charge/=gain;
116       if(charge<cal->GetThresholdAnode(ian)) continue;
117       Int_t q=(Int_t)(charge+0.5);
118       Int_t y=itb+1;
119       Int_t z=ian+1;
120       Int_t iindex=y*fNZbins+z;
121       Float_t noise=cal->GetNoiseAfterElectronics(ian)*2.2; // applies zero suppression using the measured noise of each anode. Threshold values from ALICE-INT-1999-28 V10
122       if (z<=fNAnodes){
123         if(q>noise){
124           bins0.push_back(iindex);
125           bins0.push_back(q);
126           bins0.push_back(0);
127           bins0.push_back(i);
128           map0[iindex/32]|=(1<<(iindex%32));
129         }
130       }
131       else{
132         z-=fNAnodes;
133         if(q>noise){
134           iindex=y*fNZbins+z;
135           bins1.push_back(iindex);
136           bins1.push_back(q);
137           bins1.push_back(0);
138           bins1.push_back(i);
139           map1[iindex/32]|=(1<<(iindex%32));
140         }
141       }
142     }
143   }
144   FindClustersSDD(bins0, bins1, map0, map1, digits);
145 }
146
147 //______________________________________________________________________
148 void AliITSClusterFinderSDDfast::FindClustersSDD(std::vector<int>& bins0, std::vector<int>& bins1, const Int_t map0[], const Int_t map1[], TClonesArray *digits, TClonesArray *clusters, Int_t jitter){
149
150   static AliITSRecoParam *repa = NULL;
151   if(!repa){
152     repa = (AliITSRecoParam*) AliITSReconstructor::GetRecoParam();
153     if(!repa){
154       repa = AliITSRecoParam::GetHighFluxParam();
155       AliWarning("Using default AliITSRecoParam class");
156     }
157   }
158   const TGeoHMatrix *mT2L=AliITSgeomTGeo::GetTracking2LocalMatrix(fModule);
159   AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(fModule);
160   if(cal==0){
161     AliError(Form("Calibration object not present for SDD module %d\n",fModule));
162     return;
163   }
164
165   TClonesArray &cl=*clusters;
166   Int_t nrp=0;
167   for (Int_t s=0; s<2; s++){
168     Int_t *bins;
169     unsigned int binssize;
170     const Int_t* map;
171     if(s==0){
172       binssize=bins0.size();
173       bins = &bins0[0];
174       map=map0;
175     }
176     if(s==1){
177       binssize=bins1.size();
178       bins=&bins1[0];
179       map=map1;
180     }
181
182     const Int_t rresto=fNZbins-1;
183     Int_t cid=0;
184     for(std::vector<int>::size_type i=0;i<binssize;i+=4){
185       if(!bins[i+2])
186         bins[i+2]=++cid;
187       Int_t me=bins[i];
188       Int_t resto=me%fNZbins;
189       if(resto==rresto){
190         Int_t idxs[1]={me+fNZbins};
191         for(Int_t k=0;k<1;++k)
192           if(map[idxs[k]/32]&(1<<(idxs[k]%32)))
193             for(std::vector<int>::size_type j=i+4;j<binssize;j+=4)
194               if(bins[j]==idxs[k]){
195                 bins[j+2]=bins[i+2];
196                 break;
197               }
198       }
199       else{
200         Int_t idxs[2]={me+1,me+fNZbins};
201         for(int k=0;k<2;++k)
202           if(map[idxs[k]/32]&(1<<(idxs[k]%32)))
203             for(std::vector<int>::size_type j=i+4;j<binssize;j+=4)
204               if(bins[j]==idxs[k]){
205                 bins[j+2]=bins[i+2];
206                 break;
207               }
208       }
209     }
210     for(std::vector<int>::size_type i=0;i<binssize;i+=4){
211       Int_t me=bins[i];
212       Int_t resto=me%fNZbins;
213       if(resto==fNZbins-1){
214         Int_t idxs[1]={me+fNZbins};
215         Int_t myid=bins[i+2];
216         for(Int_t l=0;l<1;++l){
217           if(map[idxs[l]/32]&(1<<(idxs[l]%32)))
218             for(std::vector<int>::size_type j=i+4;j<binssize;j+=4){
219               if(bins[j]==idxs[l]){
220                 Int_t hisid=bins[j+2];
221                 if(myid!=hisid){
222                   for(std::vector<int>::size_type k=2;k<binssize;k+=4)
223                     if(bins[k]==hisid)
224                       bins[k]=myid;
225                 }
226                 break;
227               }
228             }
229         }
230       }
231       else{
232         Int_t idxs[2]={me+1,me+fNZbins};
233         Int_t myid=bins[i+2];
234         for(Int_t l=0;l<2;++l){
235           if(map[idxs[l]/32]&(1<<(idxs[l]%32)))
236             for(std::vector<int>::size_type j=i+4;j<binssize;j+=4){
237               if(bins[j]==idxs[l]){
238                 Int_t hisid=bins[j+2];
239                 if(myid!=hisid){
240                   for(std::vector<int>::size_type k=2;k<binssize;k+=4)
241                     if(bins[k]==hisid)
242                       bins[k]=myid;
243                 }
244                 break;
245               }
246             }
247         }
248       }
249     }
250
251     Int_t recp[cid][12];
252     for(Int_t i=0;i<cid;++i)
253       for(Int_t j=0;j<12;++j)
254         recp[i][j]=0;
255
256     Int_t kplab[cid][10];
257     for(Int_t i=0;i<cid;++i)
258       for(Int_t j=0;j<10;++j)
259         kplab[i][j]=-2;
260
261     for(std::vector<int>::size_type i=0;i<binssize;i+=4){
262       Int_t q=bins[i+1];
263       Int_t me=bins[i+2]-1;
264       Int_t z=bins[i]%fNZbins;
265       Int_t x=bins[i]/fNZbins;
266       recp[me][0]+=q;                   //sumq
267       recp[me][1]+=z*q;                 //sumz
268       recp[me][2]+=x*q;                 //sumx
269
270 #ifdef CSBASEDERROR
271       recp[me][3]+=z*z*q;               //sigmaZ2
272       recp[me][4]+=x*x*q;               //sigmaX2
273 #endif
274
275       if(recp[me][5]==0){
276         recp[me][6]=z;
277         recp[me][7]=z;
278         recp[me][8]=x;
279         recp[me][9]=x;
280         recp[me][10]=q;
281         recp[me][11]=bins[i];
282       }
283       else{
284         if(recp[me][6]<z) recp[me][6]=z;
285         if(recp[me][7]>z) recp[me][7]=z;
286         if(recp[me][8]<x) recp[me][8]=x;
287         if(recp[me][9]>x) recp[me][9]=x;
288         if(recp[me][10]<q){
289           recp[me][10]=q;
290           recp[me][11]=bins[i];
291         }
292       }
293
294       if(digits){
295         Int_t kplab2[10];
296         for(Int_t ilab=0;ilab<10;++ilab)
297           kplab2[ilab]=kplab[me][ilab];
298         AliITSdigitSDD* d=(AliITSdigitSDD*)digits->UncheckedAt(bins[i+3]);
299         for (Int_t itrack=0;itrack<10;itrack++){
300           Int_t track = (d->GetTracks())[itrack];
301           if (track>=0) {
302             AddLabel(kplab2, track);
303           }
304         }
305         for(Int_t ilab=0;ilab<10;++ilab)
306           kplab[me][ilab]=kplab2[ilab];
307       }
308       ++recp[me][5];                    //nPiInClu
309     }
310
311     for(Int_t i=0;i<cid;++i){
312       if(recp[i][5]==0) continue;
313       if(recp[i][5]==1) continue;
314
315       Float_t q=recp[i][0];
316
317       Int_t clSizAnode=recp[i][6]-recp[i][7]+1;
318       Int_t clSizTb=recp[i][8]-recp[i][9]+1;
319       if(repa->GetUseSDDClusterSizeSelection()){
320         if(clSizTb==1) continue; // cut common mode noise spikes
321         if(clSizAnode>5) continue; // cut common mode noise spikes
322         if(clSizTb>10) continue; // cut clusters on noisy anodes
323         if(cal->IsAMAt20MHz() && clSizTb>8) continue; // cut clusters on noisy anodes
324       }
325
326       Float_t zz=(Float_t)recp[i][1]/q;
327       Float_t xx=(Float_t)recp[i][2]/q;
328
329       AliITSresponseSDD* rsdd = fDetTypeRec->GetResponseSDD();
330
331       Float_t zAnode=zz-0.5;  // to have anode in range 0.-255. and centered on the mid of the pitch
332       Float_t timebin=xx-0.5;  // to have time bin in range 0.-255. amd centered on the mid of the bin
333
334       if(s==1) zAnode+=fNAnodes;  // right side has anodes from 256. to 511.
335
336       Float_t zdet=GetSeg()->GetLocalZFromAnode(zAnode);
337       Float_t driftTimeUncorr=GetSeg()->GetDriftTimeFromTb(timebin)+jitter*rsdd->GetCarlosRXClockPeriod();
338       Float_t driftTime=driftTimeUncorr-rsdd->GetTimeZero(fModule);
339
340       if(driftTime<fMaxDrTimeForTightCut && recp[i][10]<fCutOnPeakTight) continue;
341
342       Float_t driftSpeed=cal->GetDriftSpeedAtAnode(zAnode) + rsdd->GetDeltaVDrift(fModule,zAnode>255);
343       Float_t driftPathMicron=driftTime*driftSpeed;
344       const Double_t kMicronTocm=1.0e-4;
345       Float_t xdet=(driftPathMicron-GetSeg()->Dx())*kMicronTocm; // xdet is negative
346       if(s==0) xdet=-xdet; // left side has positive local x
347
348       if(repa->GetUseSDDCorrectionMaps()){
349         Float_t corrx=0, corrz=0;
350         cal->GetCorrections(zdet,xdet,corrz,corrx,GetSeg());
351         zdet+=corrz;
352         xdet+=corrx;
353       }
354
355       Double_t loc[3]={xdet,0.,zdet},trk[3]={0.,0.,0.};
356       mT2L->MasterToLocal(loc,trk);
357       xx=trk[1];
358       zz=trk[2];
359
360       q+=(driftTime*rsdd->GetADCvsDriftTime(fModule)); // correction for zero supp.
361       q/=rsdd->GetADCtokeV(fModule);
362       if(cal-> IsAMAt20MHz()) q*=2.; // account for 1/2 sampling freq.
363       if(q<repa->GetMinClusterChargeSDD()) continue; // remove noise clusters
364
365 #ifdef  CSBASEDERROR
366       Float_t hit[6]={xx,zz,recp[i][3],recp[i][4],q,0.};
367 #else
368       Float_t hit[6]={xx, zz, 0.0030*0.0030, 0.0020*0.0020, q, 0.};
369 #endif
370
371       Int_t  info[3]={clSizTb, clSizAnode, fNlayer[fModule]};
372
373       Int_t kplab2[10];
374       if(digits){
375         for(Int_t ilab=0;ilab<10;++ilab)
376           if(kplab[i][ilab]!=-2)
377             kplab2[ilab]=kplab[i][ilab];
378           else
379             kplab2[ilab]=-2;
380       }
381       else{
382         if(fRawID2ClusID) kplab2[0]=fNClusters+1; // RS: store clID+1 as a reference to the cluster
383         for(Int_t ilab=1;ilab<10;++ilab)
384           kplab2[ilab]=-2;
385       }
386       if(digits) CheckLabels2(kplab2);
387       kplab2[3]=fNdet[fModule];
388       AliITSRecPoint cc(kplab2,hit,info);
389       cc.SetType(101);
390       cc.SetDriftTime(driftTimeUncorr);
391       cc.SetDriftSide(s);
392       cc.SetChargeRatio(recp[i][10]);
393       if(clusters) new (cl[nrp]) AliITSRecPoint(cc);
394       else {
395         fDetTypeRec->AddRecPoint(cc);
396       }
397       ++nrp;
398       ++fNClusters;
399     }
400   }
401   AliDebug(2,Form("Clusters found on SDD module %d (unfolding %d) = %d\n",fModule,repa->GetUseUnfoldingInClusterFinderSDD(),nrp));
402 }
403
404 //______________________________________________________________________
405 void AliITSClusterFinderSDDfast::RawdataToClusters(AliRawReader* rawReader){
406   //------------------------------------------------------------
407   // This function creates ITS clusters from raw data
408   //------------------------------------------------------------
409   fNClusters = 0; //RS
410   AliITSRawStream* inputSDD=AliITSRawStreamSDD::CreateRawStreamSDD(rawReader);
411   AliDebug(1,Form("%s is used",inputSDD->ClassName()));
412
413   AliITSDDLModuleMapSDD *ddlmap=(AliITSDDLModuleMapSDD*)fDetTypeRec->GetDDLModuleMapSDD();
414   inputSDD->SetDDLModuleMap(ddlmap);
415   for(Int_t iddl=0; iddl<AliITSDDLModuleMapSDD::GetNDDLs(); iddl++){
416     for(Int_t icar=0; icar<AliITSDDLModuleMapSDD::GetNModPerDDL();icar++){
417       Int_t iMod=ddlmap->GetModuleNumber(iddl,icar);
418       if(iMod==-1) continue;
419       AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(iMod);
420       if(cal==0){
421         AliError(Form("Calibration object not present for SDD module %d\n",iMod));
422         continue;
423       }
424       Bool_t isZeroSupp=cal->GetZeroSupp();
425       if(isZeroSupp){
426         for(Int_t iSid=0; iSid<2; iSid++) inputSDD->SetZeroSuppLowThreshold(iMod-240,iSid,cal->GetZSLowThreshold(iSid));
427       }else{
428         for(Int_t iSid=0; iSid<2; iSid++) inputSDD->SetZeroSuppLowThreshold(iMod-240,iSid,0);
429       }
430     }
431   }
432   FindClustersSDD(inputSDD);
433   delete inputSDD;
434 }
435
436 void AliITSClusterFinderSDDfast::FindClustersSDD(AliITSRawStream* input){
437
438   AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
439   Int_t nClustersSDD=0;
440
441   const Int_t kMapDim=fNZbins*fNXbins/32;
442   Int_t mapsDDL[kHybridsPerDDL][kMapDim];
443   for(Int_t i=0;i<kHybridsPerDDL;++i)
444     for(Int_t j=0;j<kMapDim;++j)
445       mapsDDL[i][j]=0;
446
447   Int_t vectModId[kModulesPerDDL];
448   for(Int_t iMod=0; iMod<kModulesPerDDL; iMod++) vectModId[iMod]=-1;
449   // read raw data input stream
450   int countRW = 0; //RS
451   if (fRawID2ClusID) fRawID2ClusID->Reset(); //RS if array was provided, we shall store the rawID -> ClusterID
452
453   while (input->Next()) {
454     Int_t iModule = input->GetModuleID();
455     if(iModule<0){
456       AliWarning(Form("Invalid SDD module number %d\n", iModule));
457       continue;
458     }
459     Int_t iCarlos=input->GetCarlosId();
460     Int_t iSide=input->GetChannel();
461     Int_t iHybrid=iCarlos*2+iSide;
462
463     if(input->IsCompletedModule()){
464       vectModId[iCarlos]=iModule; // store the module number
465     }
466     else if(input->IsCompletedDDL()){
467       // when all data from a DDL was read, search for clusters
468       Int_t jitter=input->GetJitter();
469       for(Int_t iMod=0;iMod<kModulesPerDDL;iMod++){
470         if(vectModId[iMod]>=0){
471           fModule = vectModId[iMod];
472           TClonesArray* clusters=rpc->UncheckedGetClusters(fModule);
473           std::vector<int> bins0;
474           std::vector<int> bins1;
475           bins0=fDDLBins[iMod*2];
476           bins1=fDDLBins[iMod*2+1];
477           Int_t map0[kMapDim];
478           Int_t map1[kMapDim];
479           for(Int_t j=0;j<kMapDim;++j){
480             map0[j]=map1[j]=0;
481           }
482           for(Int_t i=0;i<kMapDim;++i){
483             map0[i]=mapsDDL[iMod*2][i];
484             map1[i]=mapsDDL[iMod*2+1][i];
485           }
486
487           FindClustersSDD(bins0, bins1, map0, map1, NULL, clusters,jitter);
488
489           Int_t nClusters = clusters->GetEntriesFast();
490           nClustersSDD += nClusters;
491           vectModId[iMod]=-1;
492         }
493         for (Int_t s=0; s<2; s++){
494           Int_t indexHyb=iMod*2+s;
495           for(std::vector<int>::size_type i=0;i<fDDLBins[indexHyb].size();++i)
496             fDDLBins[indexHyb].erase(fDDLBins[indexHyb].begin(),fDDLBins[indexHyb].end());
497           for(Int_t j=0;j<kMapDim;++j)
498             mapsDDL[indexHyb][j]=0;
499         }
500       }
501     }
502     else{ // fill the current digit into the bins array
503       if(iHybrid<0 || iHybrid>=kHybridsPerDDL){
504         AliWarning(Form("Invalid SDD hybrid number %d on module %d\n", iHybrid,iModule));
505         continue;
506       }
507       AliITSCalibrationSDD* cal=(AliITSCalibrationSDD*)GetResp(iModule);
508       if(cal==0){
509         AliError(Form("Calibration object not present for SDD module %d\n",iModule));
510         continue;
511       }
512       Float_t charge=input->GetSignal();
513       Int_t chan=input->GetCoord1()+fNAnodes*iSide;
514       Float_t gain=cal->GetChannelGain(chan)/fDetTypeRec->GetAverageGainSDD();;
515       Float_t baseline=cal->GetBaseline(chan);
516       if(charge>baseline) charge-=baseline;
517       else charge=0;
518       if(gain>0.){ // Bad channels have gain=0
519         charge/=gain;
520         if(charge>=cal->GetThresholdAnode(chan)){
521           Int_t q=(Int_t)(charge+0.5);
522           Int_t iz = input->GetCoord1();
523           Int_t itb = input->GetCoord2();
524           Float_t noise=cal->GetNoiseAfterElectronics(iz)*2.2;  // applies zero suppression using the measured noise of each anode. Threshold values from ALICE-INT-1999-28 V10
525           Int_t index=(itb+1)*fNZbins+(iz+1);
526           if((itb<fNTimeBins) && (iz<fNAnodes)){
527             if(q<noise) continue;
528             fDDLBins[iHybrid].push_back(index);
529             fDDLBins[iHybrid].push_back(q);
530             fDDLBins[iHybrid].push_back(0);
531             fDDLBins[iHybrid].push_back(countRW);
532             mapsDDL[iHybrid][index/32]|=(1<<(index%32));
533           }
534           else{
535             AliWarning(Form("Invalid SDD cell: Anode=%d   TimeBin=%d",iz,itb));
536           }
537         }
538       }
539     }
540     countRW++; //RS
541   }
542   AliDebug(1,Form("found clusters in ITS SDD: %d", nClustersSDD));
543 }