]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSClusterFinderSDDfast.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Protection against missing timestamps
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSDDfast.cxx
CommitLineData
05b25e73 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
46ClassImp(AliITSClusterFinderSDDfast)
47
48AliITSClusterFinderSDDfast::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//______________________________________________________________________
70AliITSClusterFinderSDDfast::~AliITSClusterFinderSDDfast()
71{
72 //Default destructor
73}
74
75//______________________________________________________________________
76void AliITSClusterFinderSDDfast::FindRawClusters(Int_t mod){
77
78 //Find clusters
79 SetModule(mod);
80 FindClustersSDD(fDigits);
81
82}
83
84//______________________________________________________________________
85void 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//______________________________________________________________________
148void 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;
1c558e0c 184 for(std::vector<int>::size_type i=0;i<binssize;i+=4){
05b25e73 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)))
1c558e0c 193 for(std::vector<int>::size_type j=i+4;j<binssize;j+=4)
05b25e73 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)))
1c558e0c 203 for(std::vector<int>::size_type j=i+4;j<binssize;j+=4)
05b25e73 204 if(bins[j]==idxs[k]){
205 bins[j+2]=bins[i+2];
206 break;
207 }
208 }
209 }
1c558e0c 210 for(std::vector<int>::size_type i=0;i<binssize;i+=4){
05b25e73 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)))
1c558e0c 218 for(std::vector<int>::size_type j=i+4;j<binssize;j+=4){
05b25e73 219 if(bins[j]==idxs[l]){
220 Int_t hisid=bins[j+2];
221 if(myid!=hisid){
1c558e0c 222 for(std::vector<int>::size_type k=2;k<binssize;k+=4)
05b25e73 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)))
1c558e0c 236 for(std::vector<int>::size_type j=i+4;j<binssize;j+=4){
05b25e73 237 if(bins[j]==idxs[l]){
238 Int_t hisid=bins[j+2];
239 if(myid!=hisid){
1c558e0c 240 for(std::vector<int>::size_type k=2;k<binssize;k+=4)
05b25e73 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
1c558e0c 261 for(std::vector<int>::size_type i=0;i<binssize;i+=4){
05b25e73 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;
ef658336 341 if(recp[i][10]<fCutOnPeakLoose) continue;
05b25e73 342
343 Float_t driftSpeed=cal->GetDriftSpeedAtAnode(zAnode) + rsdd->GetDeltaVDrift(fModule,zAnode>255);
344 Float_t driftPathMicron=driftTime*driftSpeed;
345 const Double_t kMicronTocm=1.0e-4;
346 Float_t xdet=(driftPathMicron-GetSeg()->Dx())*kMicronTocm; // xdet is negative
347 if(s==0) xdet=-xdet; // left side has positive local x
348
349 if(repa->GetUseSDDCorrectionMaps()){
350 Float_t corrx=0, corrz=0;
351 cal->GetCorrections(zdet,xdet,corrz,corrx,GetSeg());
352 zdet+=corrz;
353 xdet+=corrx;
354 }
355
356 Double_t loc[3]={xdet,0.,zdet},trk[3]={0.,0.,0.};
357 mT2L->MasterToLocal(loc,trk);
358 xx=trk[1];
359 zz=trk[2];
360
361 q+=(driftTime*rsdd->GetADCvsDriftTime(fModule)); // correction for zero supp.
362 q/=rsdd->GetADCtokeV(fModule);
363 if(cal-> IsAMAt20MHz()) q*=2.; // account for 1/2 sampling freq.
364 if(q<repa->GetMinClusterChargeSDD()) continue; // remove noise clusters
365
366#ifdef CSBASEDERROR
367 Float_t hit[6]={xx,zz,recp[i][3],recp[i][4],q,0.};
368#else
369 Float_t hit[6]={xx, zz, 0.0030*0.0030, 0.0020*0.0020, q, 0.};
370#endif
371
372 Int_t info[3]={clSizTb, clSizAnode, fNlayer[fModule]};
373
374 Int_t kplab2[10];
375 if(digits){
376 for(Int_t ilab=0;ilab<10;++ilab)
377 if(kplab[i][ilab]!=-2)
378 kplab2[ilab]=kplab[i][ilab];
379 else
380 kplab2[ilab]=-2;
381 }
382 else{
383 if(fRawID2ClusID) kplab2[0]=fNClusters+1; // RS: store clID+1 as a reference to the cluster
384 for(Int_t ilab=1;ilab<10;++ilab)
385 kplab2[ilab]=-2;
386 }
387 if(digits) CheckLabels2(kplab2);
388 kplab2[3]=fNdet[fModule];
389 AliITSRecPoint cc(kplab2,hit,info);
390 cc.SetType(101);
391 cc.SetDriftTime(driftTimeUncorr);
392 cc.SetDriftSide(s);
393 cc.SetChargeRatio(recp[i][10]);
394 if(clusters) new (cl[nrp]) AliITSRecPoint(cc);
395 else {
396 fDetTypeRec->AddRecPoint(cc);
397 }
398 ++nrp;
399 ++fNClusters;
400 }
401 }
402 AliDebug(2,Form("Clusters found on SDD module %d (unfolding %d) = %d\n",fModule,repa->GetUseUnfoldingInClusterFinderSDD(),nrp));
403}
404
405//______________________________________________________________________
406void AliITSClusterFinderSDDfast::RawdataToClusters(AliRawReader* rawReader){
407 //------------------------------------------------------------
408 // This function creates ITS clusters from raw data
409 //------------------------------------------------------------
410 fNClusters = 0; //RS
411 AliITSRawStream* inputSDD=AliITSRawStreamSDD::CreateRawStreamSDD(rawReader);
412 AliDebug(1,Form("%s is used",inputSDD->ClassName()));
413
414 AliITSDDLModuleMapSDD *ddlmap=(AliITSDDLModuleMapSDD*)fDetTypeRec->GetDDLModuleMapSDD();
415 inputSDD->SetDDLModuleMap(ddlmap);
416 for(Int_t iddl=0; iddl<AliITSDDLModuleMapSDD::GetNDDLs(); iddl++){
417 for(Int_t icar=0; icar<AliITSDDLModuleMapSDD::GetNModPerDDL();icar++){
418 Int_t iMod=ddlmap->GetModuleNumber(iddl,icar);
419 if(iMod==-1) continue;
420 AliITSCalibrationSDD* cal = (AliITSCalibrationSDD*)GetResp(iMod);
421 if(cal==0){
422 AliError(Form("Calibration object not present for SDD module %d\n",iMod));
423 continue;
424 }
425 Bool_t isZeroSupp=cal->GetZeroSupp();
426 if(isZeroSupp){
427 for(Int_t iSid=0; iSid<2; iSid++) inputSDD->SetZeroSuppLowThreshold(iMod-240,iSid,cal->GetZSLowThreshold(iSid));
428 }else{
429 for(Int_t iSid=0; iSid<2; iSid++) inputSDD->SetZeroSuppLowThreshold(iMod-240,iSid,0);
430 }
431 }
432 }
433 FindClustersSDD(inputSDD);
434 delete inputSDD;
435}
436
437void AliITSClusterFinderSDDfast::FindClustersSDD(AliITSRawStream* input){
438
439 AliITSRecPointContainer* rpc = AliITSRecPointContainer::Instance();
440 Int_t nClustersSDD=0;
441
442 const Int_t kMapDim=fNZbins*fNXbins/32;
443 Int_t mapsDDL[kHybridsPerDDL][kMapDim];
444 for(Int_t i=0;i<kHybridsPerDDL;++i)
445 for(Int_t j=0;j<kMapDim;++j)
446 mapsDDL[i][j]=0;
447
448 Int_t vectModId[kModulesPerDDL];
449 for(Int_t iMod=0; iMod<kModulesPerDDL; iMod++) vectModId[iMod]=-1;
450 // read raw data input stream
451 int countRW = 0; //RS
452 if (fRawID2ClusID) fRawID2ClusID->Reset(); //RS if array was provided, we shall store the rawID -> ClusterID
453
454 while (input->Next()) {
455 Int_t iModule = input->GetModuleID();
456 if(iModule<0){
457 AliWarning(Form("Invalid SDD module number %d\n", iModule));
458 continue;
459 }
460 Int_t iCarlos=input->GetCarlosId();
461 Int_t iSide=input->GetChannel();
462 Int_t iHybrid=iCarlos*2+iSide;
463
464 if(input->IsCompletedModule()){
465 vectModId[iCarlos]=iModule; // store the module number
466 }
467 else if(input->IsCompletedDDL()){
468 // when all data from a DDL was read, search for clusters
469 Int_t jitter=input->GetJitter();
470 for(Int_t iMod=0;iMod<kModulesPerDDL;iMod++){
471 if(vectModId[iMod]>=0){
472 fModule = vectModId[iMod];
473 TClonesArray* clusters=rpc->UncheckedGetClusters(fModule);
1c558e0c 474 std::vector<int> bins0;
475 std::vector<int> bins1;
05b25e73 476 bins0=fDDLBins[iMod*2];
477 bins1=fDDLBins[iMod*2+1];
478 Int_t map0[kMapDim];
479 Int_t map1[kMapDim];
480 for(Int_t j=0;j<kMapDim;++j){
481 map0[j]=map1[j]=0;
482 }
483 for(Int_t i=0;i<kMapDim;++i){
484 map0[i]=mapsDDL[iMod*2][i];
485 map1[i]=mapsDDL[iMod*2+1][i];
486 }
487
488 FindClustersSDD(bins0, bins1, map0, map1, NULL, clusters,jitter);
489
490 Int_t nClusters = clusters->GetEntriesFast();
491 nClustersSDD += nClusters;
492 vectModId[iMod]=-1;
493 }
494 for (Int_t s=0; s<2; s++){
495 Int_t indexHyb=iMod*2+s;
1c558e0c 496 for(std::vector<int>::size_type i=0;i<fDDLBins[indexHyb].size();++i)
05b25e73 497 fDDLBins[indexHyb].erase(fDDLBins[indexHyb].begin(),fDDLBins[indexHyb].end());
498 for(Int_t j=0;j<kMapDim;++j)
499 mapsDDL[indexHyb][j]=0;
500 }
501 }
502 }
503 else{ // fill the current digit into the bins array
504 if(iHybrid<0 || iHybrid>=kHybridsPerDDL){
505 AliWarning(Form("Invalid SDD hybrid number %d on module %d\n", iHybrid,iModule));
506 continue;
507 }
508 AliITSCalibrationSDD* cal=(AliITSCalibrationSDD*)GetResp(iModule);
509 if(cal==0){
510 AliError(Form("Calibration object not present for SDD module %d\n",iModule));
511 continue;
512 }
513 Float_t charge=input->GetSignal();
514 Int_t chan=input->GetCoord1()+fNAnodes*iSide;
515 Float_t gain=cal->GetChannelGain(chan)/fDetTypeRec->GetAverageGainSDD();;
516 Float_t baseline=cal->GetBaseline(chan);
517 if(charge>baseline) charge-=baseline;
518 else charge=0;
519 if(gain>0.){ // Bad channels have gain=0
520 charge/=gain;
521 if(charge>=cal->GetThresholdAnode(chan)){
522 Int_t q=(Int_t)(charge+0.5);
523 Int_t iz = input->GetCoord1();
524 Int_t itb = input->GetCoord2();
ef658336 525 Int_t ian=iz;
526 if(iSide==1)
527 ian+=256;
528 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
05b25e73 529 Int_t index=(itb+1)*fNZbins+(iz+1);
530 if((itb<fNTimeBins) && (iz<fNAnodes)){
531 if(q<noise) continue;
ef658336 532 fDDLBins[iHybrid].push_back(index);
533 fDDLBins[iHybrid].push_back(q);
534 fDDLBins[iHybrid].push_back(0);
535 fDDLBins[iHybrid].push_back(countRW);
536 mapsDDL[iHybrid][index/32]|=(1<<(index%32));
05b25e73 537 }
538 else{
539 AliWarning(Form("Invalid SDD cell: Anode=%d TimeBin=%d",iz,itb));
540 }
541 }
542 }
543 }
544 countRW++; //RS
545 }
546 AliDebug(1,Form("found clusters in ITS SDD: %d", nClustersSDD));
547}