V2 clusterer moved to the standard framework. V2 clusters and recpoints are still...
[u/mrichter/AliRoot.git] / ITS / AliITSClusterFinderSDD.cxx
CommitLineData
b0f5e3fc 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 **************************************************************************/
d86f531c 15/*
16 $Id$
17 $Log$
04366a57 18 Revision 1.37 2004/06/10 21:00:24 nilsen
19 Modifications associated with remerging the Ba/Sa and Dubna pixel simulations,
20 some cleaning of general code (including coding convensions), and adding some
21 protections associated with SetDefaults/SetDefaultSimulations which should help
22 with the Test beam simulations. Details below. The default SPD simulation for
23 the general ITS runs/geometry is still the Ba/Sa, but for the Test beam
24 geometries this has been changed to the merged versions.
25 File: AliITS.cxx Modified
26 File: AliITS.h Modified
27 In lined many one-two line functions. Added some protection to
28 SetDefaults(), SetDefaultSimulation(), and SetDefaultClusterFinders(),
29 such that they should now even work when only one detector type has
30 been defined (as it should be for the test beams...). Some mostly
31 cosmetic issues associated with getting branch names for digits. And
32 Generally some cleaning up of the code.
33 File: AliITSClusterFinder.cxx Modified
34 File: AliITSClusterFinder.h Modified
35 Did some additional consolidation of data into the base class, added
36 TClonesArray *fClusters, a fDebug, and fModule variables. Otherwise
37 some cosmetic and coding conversion changes.
38 File: AliITSClusterFinderSDD.cxx Modified
39 File: AliITSClusterFinderSDD.h Modified
40 Changes to be consistent with the modified base class, and cosmetic
41 and coding conversion changes.
42 File: AliITSClusterFinderSPD.cxx Modified
43 File: AliITSClusterFinderSPD.h Modified
44 Changes to be consistent with the modified base class, and cosmetic
45 and coding conversion changes.
46 File: AliITSClusterFinderSPDdubna.h Removed
47 File: AliITSClusterFinderSPDdubna.cxx Removed
48 Since we have ClusterFinderSPD and V2 and this version isn't being
49 maintained, it is being retired.
50 File: AliITSClusterFinderSSD.cxx Modified
51 File: AliITSClusterFinderSSD.h Modified
52 Changes to be consistent with the modified base class, and cosmetic
53 and coding conversion changes.
54 File: AliITSDetType.cxx Modified
55 File: AliITSDetType.h Modified
56 Added a new class variable to indicate what the detector type is
57 AliITSDetector fDetType; values of kSPD, kSDD, kSSD, .... Otherwise
58 cosmetic and Coding convention changes.
59 File: AliITSLoader.cxx Modified
60 File: AliITSLoader.h Modified
61 Some changes which are not complete. The idea is to be able to get,
62 simply via one call, a specific hit, Sdigit, digit, RecPoint,...
63 without all of the usual over head of initializing TClonesArrays setting
64 branch addresses and the like. Work is far form ready.
65 File: AliITSdcsSSD.cxx Modified
66 Some nearly cosmetic changes necessary due to changes to response and
67 segmentation class'.
68 File: AliITSgeom.h Modified
69 In the definition of AliITSDetector type, added kND=-1, no detector
70 defined. Expect to use it later(?).
71 File: AliITSresponse.h Modified
72 Basically cosmetic. Mostly changing Float_t to Double_t.
73 File: AliITSresponseSDD.cxx Modified
74 File: AliITSresponseSDD.h Modified
75 Basically the cosmetic and Float_t to Double_t
76 File: AliITSresponseSPD.cxx Modified
77 File: AliITSresponseSPD.h Modified
78 Mostly Float_t to Double_t and added in the IsPixelDead function for
79 the dubna version (otherwise the merging had been done).
80 File: AliITSresponseSPDdubna.h Removed
81 File: AliITSresponseSPDdubna.cxx Removed
82 We should be able to remove this class now. AliITSresponseSPD is now
83 used for both the Bari-Salerno and the dubna models.
84 File: AliITSresponseSSD.cxx Modified
85 File: AliITSresponseSSD.h Modified
86 Float_t to Double_t changes.
87 File: AliITSsegmentation.h Modified
88 Made LocaltoDet return a Bool_t. Now if the x,z location is outside
89 of the volume, it returns kFALSE. see below.
90 File: AliITSsegmentationSDD.cxx Modified
91 File: AliITSsegmentationSDD.h Modified
92 Made LocaltoDet return a Bool_t. Now if the x,z location is outside
93 of the volume, it returns kFALSE.
94 File: AliITSsegmentationSPD.cxx Modified
95 File: AliITSsegmentationSPD.h Modified
96 Made LocaltoDet return a Bool_t. Now if the x,z location is outside
97 of the volume, it returns kFALSE.
98 File: AliITSsegmentationSSD.cxx Modified
99 File: AliITSsegmentationSSD.h Modified
100 Made LocaltoDet return a Bool_t. Now if the x,z location is outside
101 of the volume, it returns kFALSE. see below.
102 File: AliITSsimulation.cxx Modified
103 File: AliITSsimulation.h Modified
104 Added fDebug variable, new Constructor for use below. Cosmetic and
105 coding convention changes
106 File: AliITSsimulationSDD.cxx Modified
107 File: AliITSsimulationSDD.h Modified
108 Added new Constructor, removed redundant variables and Cosmetic and
109 coding convention changes.
110 File: AliITSsimulationSPD.cxx Modified
111 File: AliITSsimulationSPD.h Modified
112 Removed some dead code, made changes as needed by the changes above
113 (response and segmentation classes...). a few cosmetic and coding
114 convention changes.
115 File: AliITSsimulationSPDdubna.cxx Modified
116 File: AliITSsimulationSPDdubna.h Modified
117 New merged version, implemented new and old coupling with switch,
118 coding convention and similar changes. (found 1 bugs, missing
119 ! in front of if(mod-LineSegmentL(....,).
120 File: AliITSsimulationSSD.cxx Modified
121 File: AliITSsimulationSSD.h Modified
122 removed redundant variables with base class. Fixed for coding
123 convention and other cosmetic changes.
124 File: AliITSvSDD03.cxx Modified
125 File: AliITSvSPD02.cxx Modified
126 File: AliITSvSSD03.cxx Modified
127 These two have their private versions of SetDefaults and
128 SetDefaultSimulation which have been similarly protected as in AliITS.cxx
129 File: ITSLinkDef.h Modified
130 File: libITS.pkg Modified
131 Versions which include v11 geometry and other private changes
132
aacedc3e 133 Revision 1.36 2004/01/27 16:12:03 masera
134 Coding conventions for AliITSdigitXXX classes and AliITSTrackerV1
135
e869281d 136 Revision 1.35 2003/11/10 16:33:50 masera
137 Changes to obey our coding conventions
138
41b19549 139 Revision 1.34 2003/09/11 13:48:52 masera
140 Data members of AliITSdigit classes defined as protected (They were public)
141
ecee53fc 142 Revision 1.33 2003/07/21 14:20:51 masera
143 Fix to track labes in SDD Rec-points
144
d86f531c 145 Revision 1.31.2.1 2003/07/16 13:18:04 masera
146 Proper fix to track labels associated to SDD rec-points
507cf1a7 147
d86f531c 148 Revision 1.31 2003/05/19 14:44:41 masera
149 Fix to track labels associated to SDD rec-points
bf3f2830 150
d86f531c 151 Revision 1.30 2003/03/03 16:34:35 masera
152 Corrections to comply with coding conventions
153
154 Revision 1.29 2002/10/25 18:54:22 barbera
155 Various improvements and updates from B.S.Nilsen and T. Virgili
156
157 Revision 1.28 2002/10/22 14:45:29 alibrary
158 Introducing Riostream.h
159
160 Revision 1.27 2002/10/14 14:57:00 hristov
161 Merging the VirtualMC branch to the main development branch (HEAD)
162
163 Revision 1.23.4.2 2002/10/14 13:14:07 hristov
164 Updating VirtualMC to v3-09-02
165
166 Revision 1.26 2002/09/09 17:23:28 nilsen
167 Minor changes in support of changes to AliITSdigitS?D class'.
168
169 Revision 1.25 2002/05/10 22:29:40 nilsen
170 Change my Massimo Masera in the default constructor to bring things into
171 compliance.
172
173 Revision 1.24 2002/04/24 22:02:31 nilsen
174 New SDigits and Digits routines, and related changes, (including new
175 noise values).
176
177 */
04366a57 178///////////////////////////////////////////////////////////////////////////
179// Cluster finder //
180// for Silicon //
181// Drift Detector //
182//////////////////////////////////////////////////////////////////////////
4ae5bbc4 183#include <Riostream.h>
bf3f2830 184
a1f090e0 185#include <TMath.h>
186#include <math.h>
b0f5e3fc 187
188#include "AliITSClusterFinderSDD.h"
e8189707 189#include "AliITSMapA1.h"
190#include "AliITS.h"
e869281d 191#include "AliITSdigitSDD.h"
41b19549 192#include "AliITSRawClusterSDD.h"
78a228db 193#include "AliITSRecPoint.h"
aacedc3e 194#include "AliITSsegmentationSDD.h"
5dd4cc39 195#include "AliITSresponseSDD.h"
b0f5e3fc 196#include "AliRun.h"
197
b0f5e3fc 198ClassImp(AliITSClusterFinderSDD)
199
42da2935 200//______________________________________________________________________
aacedc3e 201AliITSClusterFinderSDD::AliITSClusterFinderSDD():
202AliITSClusterFinder(),
203fNclusters(0),
204fDAnode(0.0),
205fDTime(0.0),
206fTimeCorr(0.0),
207fCutAmplitude(0),
208fMinPeak(0),
209fMinCharge(0),
210fMinNCells(0),
211fMaxNCells(0){
212 // default constructor
213}
214//______________________________________________________________________
42da2935 215AliITSClusterFinderSDD::AliITSClusterFinderSDD(AliITSsegmentation *seg,
50d05d7b 216 AliITSresponse *response,
217 TClonesArray *digits,
aacedc3e 218 TClonesArray *recp):
219AliITSClusterFinder(seg,response),
220fNclusters(0),
221fDAnode(0.0),
222fDTime(0.0),
223fTimeCorr(0.0),
224fCutAmplitude(0),
225fMinPeak(0),
226fMinCharge(0),
227fMinNCells(0),
228fMaxNCells(0){
42da2935 229 // standard constructor
78a228db 230
aacedc3e 231 SetDigits(digits);
232 SetClusters(recp);
b0f5e3fc 233 SetCutAmplitude();
234 SetDAnode();
235 SetDTime();
aacedc3e 236 SetMinPeak((Int_t)(((AliITSresponseSDD*)GetResp())->
237 GetNoiseAfterElectronics()*5));
355ccb70 238 // SetMinPeak();
78a228db 239 SetMinNCells();
240 SetMaxNCells();
241 SetTimeCorr();
a1f090e0 242 SetMinCharge();
aacedc3e 243 SetMap(new AliITSMapA1(GetSeg(),Digits(),fCutAmplitude));
e8189707 244}
42da2935 245//______________________________________________________________________
aacedc3e 246void AliITSClusterFinderSDD::SetCutAmplitude(Double_t nsigma){
42da2935 247 // set the signal threshold for cluster finder
aacedc3e 248 Double_t baseline,noise,noiseAfterEl;
42da2935 249
aacedc3e 250 GetResp()->GetNoiseParam(noise,baseline);
251 noiseAfterEl = ((AliITSresponseSDD*)GetResp())->GetNoiseAfterElectronics();
bf3f2830 252 fCutAmplitude = (Int_t)((baseline + nsigma*noiseAfterEl));
5dd4cc39 253}
42da2935 254//______________________________________________________________________
255void AliITSClusterFinderSDD::Find1DClusters(){
256 // find 1D clusters
a1f090e0 257
42da2935 258 // retrieve the parameters
aacedc3e 259 Int_t fNofMaps = GetSeg()->Npz();
260 Int_t fMaxNofSamples = GetSeg()->Npx();
b48af428 261 Int_t fNofAnodes = fNofMaps/2;
262 Int_t dummy = 0;
aacedc3e 263 Double_t fTimeStep = GetSeg()->Dpx(dummy);
264 Double_t fSddLength = GetSeg()->Dx();
265 Double_t fDriftSpeed = GetResp()->DriftSpeed();
266 Double_t anodePitch = GetSeg()->Dpz(dummy);
42da2935 267
268 // map the signal
aacedc3e 269 Map()->ClearMap();
270 Map()->SetThreshold(fCutAmplitude);
271 Map()->FillMap();
a1f090e0 272
aacedc3e 273 Double_t noise;
274 Double_t baseline;
275 GetResp()->GetNoiseParam(noise,baseline);
a1f090e0 276
42da2935 277 Int_t nofFoundClusters = 0;
278 Int_t i;
aacedc3e 279 Double_t **dfadc = new Double_t*[fNofAnodes];
280 for(i=0;i<fNofAnodes;i++) dfadc[i] = new Double_t[fMaxNofSamples];
281 Double_t fadc = 0.;
282 Double_t fadc1 = 0.;
283 Double_t fadc2 = 0.;
42da2935 284 Int_t j,k,idx,l,m;
285 for(j=0;j<2;j++) {
50d05d7b 286 for(k=0;k<fNofAnodes;k++) {
287 idx = j*fNofAnodes+k;
288 // signal (fadc) & derivative (dfadc)
289 dfadc[k][255]=0.;
290 for(l=0; l<fMaxNofSamples; l++) {
aacedc3e 291 fadc2=(Double_t)Map()->GetSignal(idx,l);
292 if(l>0) fadc1=(Double_t)Map()->GetSignal(idx,l-1);
50d05d7b 293 if(l>0) dfadc[k][l-1] = fadc2-fadc1;
294 } // samples
295 } // anodes
42da2935 296
50d05d7b 297 for(k=0;k<fNofAnodes;k++) {
aacedc3e 298 if(GetDebug(5)) cout<<"Anode: "<<k+1<<", Wing: "<<j+1<< endl;
50d05d7b 299 idx = j*fNofAnodes+k;
300 Int_t imax = 0;
301 Int_t imaxd = 0;
302 Int_t it = 0;
303 while(it <= fMaxNofSamples-3) {
304 imax = it;
305 imaxd = it;
306 // maximum of signal
aacedc3e 307 Double_t fadcmax = 0.;
308 Double_t dfadcmax = 0.;
50d05d7b 309 Int_t lthrmina = 1;
310 Int_t lthrmint = 3;
311 Int_t lthra = 1;
312 Int_t lthrt = 0;
313 for(m=0;m<20;m++) {
314 Int_t id = it+m;
315 if(id>=fMaxNofSamples) break;
aacedc3e 316 fadc=(float)Map()->GetSignal(idx,id);
50d05d7b 317 if(fadc > fadcmax) { fadcmax = fadc; imax = id;}
aacedc3e 318 if(fadc > (float)fCutAmplitude)lthrt++;
50d05d7b 319 if(dfadc[k][id] > dfadcmax) {
320 dfadcmax = dfadc[k][id];
321 imaxd = id;
322 } // end if
323 } // end for m
324 it = imaxd;
aacedc3e 325 if(Map()->TestHit(idx,imax) == kEmpty) {it++; continue;}
50d05d7b 326 // cluster charge
327 Int_t tstart = it-2;
328 if(tstart < 0) tstart = 0;
329 Bool_t ilcl = 0;
330 if(lthrt >= lthrmint && lthra >= lthrmina) ilcl = 1;
331 if(ilcl) {
332 nofFoundClusters++;
333 Int_t tstop = tstart;
aacedc3e 334 Double_t dfadcmin = 10000.;
50d05d7b 335 Int_t ij;
336 for(ij=0; ij<20; ij++) {
337 if(tstart+ij > 255) { tstop = 255; break; }
aacedc3e 338 fadc=(float)Map()->GetSignal(idx,tstart+ij);
50d05d7b 339 if((dfadc[k][tstart+ij] < dfadcmin) &&
340 (fadc > fCutAmplitude)) {
341 tstop = tstart+ij+5;
342 if(tstop > 255) tstop = 255;
343 dfadcmin = dfadc[k][it+ij];
344 } // end if
345 } // end for ij
42da2935 346
aacedc3e 347 Double_t clusterCharge = 0.;
348 Double_t clusterAnode = k+0.5;
349 Double_t clusterTime = 0.;
50d05d7b 350 Int_t clusterMult = 0;
aacedc3e 351 Double_t clusterPeakAmplitude = 0.;
50d05d7b 352 Int_t its,peakpos = -1;
aacedc3e 353 Double_t n, baseline;
354 GetResp()->GetNoiseParam(n,baseline);
50d05d7b 355 for(its=tstart; its<=tstop; its++) {
aacedc3e 356 fadc=(float)Map()->GetSignal(idx,its);
50d05d7b 357 if(fadc>baseline) fadc -= baseline;
358 else fadc = 0.;
359 clusterCharge += fadc;
360 // as a matter of fact we should take the peak
361 // pos before FFT
362 // to get the list of tracks !!!
363 if(fadc > clusterPeakAmplitude) {
364 clusterPeakAmplitude = fadc;
aacedc3e 365 //peakpos=Map()->GetHitIndex(idx,its);
50d05d7b 366 Int_t shift = (int)(fTimeCorr/fTimeStep);
367 if(its>shift && its<(fMaxNofSamples-shift))
aacedc3e 368 peakpos = Map()->GetHitIndex(idx,its+shift);
369 else peakpos = Map()->GetHitIndex(idx,its);
370 if(peakpos<0) peakpos =Map()->GetHitIndex(idx,its);
50d05d7b 371 } // end if
372 clusterTime += fadc*its;
373 if(fadc > 0) clusterMult++;
374 if(its == tstop) {
375 clusterTime /= (clusterCharge/fTimeStep); // ns
376 if(clusterTime>fTimeCorr) clusterTime -=fTimeCorr;
377 //ns
378 } // end if
379 } // end for its
42da2935 380
aacedc3e 381 Double_t clusteranodePath = (clusterAnode - fNofAnodes/2)*
382 anodePitch;
383 Double_t clusterDriftPath = clusterTime*fDriftSpeed;
50d05d7b 384 clusterDriftPath = fSddLength-clusterDriftPath;
385 if(clusterCharge <= 0.) break;
386 AliITSRawClusterSDD clust(j+1,//i
387 clusterAnode,clusterTime,//ff
388 clusterCharge, //f
389 clusterPeakAmplitude, //f
390 peakpos, //i
391 0.,0.,clusterDriftPath,//fff
392 clusteranodePath, //f
393 clusterMult, //i
394 0,0,0,0,0,0,0);//7*i
04366a57 395 fITS->AddCluster(1,&clust);
50d05d7b 396 it = tstop;
397 } // ilcl
398 it++;
399 } // while (samples)
400 } // anodes
42da2935 401 } // detectors (2)
42da2935 402
403 for(i=0;i<fNofAnodes;i++) delete[] dfadc[i];
404 delete [] dfadc;
a1f090e0 405
42da2935 406 return;
a1f090e0 407}
42da2935 408//______________________________________________________________________
409void AliITSClusterFinderSDD::Find1DClustersE(){
24a1c341 410 // find 1D clusters
42da2935 411 // retrieve the parameters
aacedc3e 412 Int_t fNofMaps = GetSeg()->Npz();
413 Int_t fMaxNofSamples = GetSeg()->Npx();
42da2935 414 Int_t fNofAnodes = fNofMaps/2;
415 Int_t dummy=0;
aacedc3e 416 Double_t fTimeStep = GetSeg()->Dpx( dummy );
417 Double_t fSddLength = GetSeg()->Dx();
418 Double_t fDriftSpeed = GetResp()->DriftSpeed();
419 Double_t anodePitch = GetSeg()->Dpz( dummy );
420 Double_t n, baseline;
421 GetResp()->GetNoiseParam( n, baseline );
42da2935 422 // map the signal
aacedc3e 423 Map()->ClearMap();
424 Map()->SetThreshold( fCutAmplitude );
425 Map()->FillMap();
50d05d7b 426
42da2935 427 Int_t nClu = 0;
50d05d7b 428 // cout << "Search cluster... "<< endl;
42da2935 429 for( Int_t j=0; j<2; j++ ){
50d05d7b 430 for( Int_t k=0; k<fNofAnodes; k++ ){
431 Int_t idx = j*fNofAnodes+k;
432 Bool_t on = kFALSE;
433 Int_t start = 0;
434 Int_t nTsteps = 0;
aacedc3e 435 Double_t fmax = 0.;
50d05d7b 436 Int_t lmax = 0;
aacedc3e 437 Double_t charge = 0.;
438 Double_t time = 0.;
439 Double_t anode = k+0.5;
50d05d7b 440 Int_t peakpos = -1;
441 for( Int_t l=0; l<fMaxNofSamples; l++ ){
aacedc3e 442 Double_t fadc = (Double_t)Map()->GetSignal( idx, l );
50d05d7b 443 if( fadc > 0.0 ){
444 if( on == kFALSE && l<fMaxNofSamples-4){
445 // star RawCluster (reset var.)
aacedc3e 446 Double_t fadc1 = (Double_t)Map()->GetSignal( idx, l+1 );
50d05d7b 447 if( fadc1 < fadc ) continue;
448 start = l;
449 fmax = 0.;
450 lmax = 0;
451 time = 0.;
452 charge = 0.;
453 on = kTRUE;
454 nTsteps = 0;
455 } // end if on...
456 nTsteps++ ;
457 if( fadc > baseline ) fadc -= baseline;
458 else fadc=0.;
459 charge += fadc;
460 time += fadc*l;
461 if( fadc > fmax ){
462 fmax = fadc;
463 lmax = l;
464 Int_t shift = (Int_t)(fTimeCorr/fTimeStep + 0.5);
465 if( l > shift && l < (fMaxNofSamples-shift) )
aacedc3e 466 peakpos = Map()->GetHitIndex( idx, l+shift );
50d05d7b 467 else
aacedc3e 468 peakpos = Map()->GetHitIndex( idx, l );
469 if( peakpos < 0) peakpos = Map()->GetHitIndex(idx,l);
50d05d7b 470 } // end if fadc
471 }else{ // end fadc>0
472 if( on == kTRUE ){
473 if( nTsteps > 2 ){
474 // min # of timesteps for a RawCluster
475 // Found a RawCluster...
476 Int_t stop = l-1;
477 time /= (charge/fTimeStep); // ns
478 // time = lmax*fTimeStep; // ns
479 if( time > fTimeCorr ) time -= fTimeCorr; // ns
aacedc3e 480 Double_t anodePath =(anode-fNofAnodes/2)*anodePitch;
481 Double_t driftPath = time*fDriftSpeed;
50d05d7b 482 driftPath = fSddLength-driftPath;
483 AliITSRawClusterSDD clust(j+1,anode,time,charge,
484 fmax, peakpos,0.,0.,
485 driftPath,anodePath,
486 nTsteps,start,stop,
487 start, stop, 1, k, k );
04366a57 488 fITS->AddCluster( 1, &clust );
aacedc3e 489 if(GetDebug(5)) clust.PrintInfo();
50d05d7b 490 nClu++;
491 } // end if nTsteps
492 on = kFALSE;
493 } // end if on==kTRUE
494 } // end if fadc>0
495 } // samples
496 } // anodes
42da2935 497 } // wings
aacedc3e 498 if(GetDebug(3)) cout << "# Rawclusters " << nClu << endl;
42da2935 499 return;
a1f090e0 500}
42da2935 501//_______________________________________________________________________
aacedc3e 502Int_t AliITSClusterFinderSDD::SearchPeak(Double_t *spect,Int_t xdim,Int_t zdim,
50d05d7b 503 Int_t *peakX, Int_t *peakZ,
aacedc3e 504 Double_t *peakAmp, Double_t minpeak ){
42da2935 505 // search peaks on a 2D cluster
506 Int_t npeak = 0; // # peaks
56fff130 507 Int_t i,j;
42da2935 508 // search peaks
509 for( Int_t z=1; z<zdim-1; z++ ){
48058160 510 for( Int_t x=1; x<xdim-2; x++ ){
aacedc3e 511 Double_t sxz = spect[x*zdim+z];
512 Double_t sxz1 = spect[(x+1)*zdim+z];
513 Double_t sxz2 = spect[(x-1)*zdim+z];
50d05d7b 514 // search a local max. in s[x,z]
515 if( sxz < minpeak || sxz1 <= 0 || sxz2 <= 0 ) continue;
516 if( sxz >= spect[(x+1)*zdim+z ] && sxz >= spect[(x-1)*zdim+z ] &&
517 sxz >= spect[x*zdim +z+1] && sxz >= spect[x*zdim +z-1] &&
518 sxz >= spect[(x+1)*zdim+z+1] && sxz >= spect[(x+1)*zdim+z-1] &&
519 sxz >= spect[(x-1)*zdim+z+1] && sxz >= spect[(x-1)*zdim+z-1] ){
520 // peak found
521 peakX[npeak] = x;
522 peakZ[npeak] = z;
523 peakAmp[npeak] = sxz;
524 npeak++;
525 } // end if ....
526 } // end for x
42da2935 527 } // end for z
528 // search groups of peaks with same amplitude.
529 Int_t *flag = new Int_t[npeak];
530 for( i=0; i<npeak; i++ ) flag[i] = 0;
531 for( i=0; i<npeak; i++ ){
50d05d7b 532 for( j=0; j<npeak; j++ ){
533 if( i==j) continue;
534 if( flag[j] > 0 ) continue;
535 if( peakAmp[i] == peakAmp[j] &&
536 TMath::Abs(peakX[i]-peakX[j])<=1 &&
537 TMath::Abs(peakZ[i]-peakZ[j])<=1 ){
538 if( flag[i] == 0) flag[i] = i+1;
539 flag[j] = flag[i];
540 } // end if ...
541 } // end for j
42da2935 542 } // end for i
50d05d7b 543 // make average of peak groups
42da2935 544 for( i=0; i<npeak; i++ ){
50d05d7b 545 Int_t nFlag = 1;
546 if( flag[i] <= 0 ) continue;
547 for( j=0; j<npeak; j++ ){
548 if( i==j ) continue;
549 if( flag[j] != flag[i] ) continue;
550 peakX[i] += peakX[j];
551 peakZ[i] += peakZ[j];
552 nFlag++;
553 npeak--;
554 for( Int_t k=j; k<npeak; k++ ){
555 peakX[k] = peakX[k+1];
556 peakZ[k] = peakZ[k+1];
557 peakAmp[k] = peakAmp[k+1];
558 flag[k] = flag[k+1];
559 } // end for k
560 j--;
561 } // end for j
562 if( nFlag > 1 ){
563 peakX[i] /= nFlag;
564 peakZ[i] /= nFlag;
565 } // end fi nFlag
42da2935 566 } // end for i
567 delete [] flag;
568 return( npeak );
a1f090e0 569}
42da2935 570//______________________________________________________________________
aacedc3e 571void AliITSClusterFinderSDD::PeakFunc( Int_t xdim, Int_t zdim, Double_t *par,
572 Double_t *spe, Double_t *integral){
24a1c341 573 // function used to fit the clusters
50d05d7b 574 // par -> parameters..
24a1c341 575 // par[0] number of peaks.
576 // for each peak i=1, ..., par[0]
50d05d7b 577 // par[i] = Ampl.
578 // par[i+1] = xpos
579 // par[i+2] = zpos
580 // par[i+3] = tau
581 // par[i+4] = sigma.
aacedc3e 582 Int_t electronics = GetResp()->Electronics(); // 1 = PASCAL, 2 = OLA
24a1c341 583 const Int_t knParam = 5;
584 Int_t npeak = (Int_t)par[0];
42da2935 585
aacedc3e 586 memset( spe, 0, sizeof( Double_t )*zdim*xdim );
42da2935 587
24a1c341 588 Int_t k = 1;
42da2935 589 for( Int_t i=0; i<npeak; i++ ){
24a1c341 590 if( integral != 0 ) integral[i] = 0.;
aacedc3e 591 Double_t sigmaA2 = par[k+4]*par[k+4]*2.;
592 Double_t t2 = par[k+3]; // PASCAL
bf3f2830 593 if( electronics == 2 ) { t2 *= t2; t2 *= 2; } // OLA
42da2935 594 for( Int_t z=0; z<zdim; z++ ){
595 for( Int_t x=0; x<xdim; x++ ){
aacedc3e 596 Double_t z2 = (z-par[k+2])*(z-par[k+2])/sigmaA2;
597 Double_t x2 = 0.;
598 Double_t signal = 0.;
42da2935 599 if( electronics == 1 ){ // PASCAL
bf3f2830 600 x2 = (x-par[k+1]+t2)/t2;
42da2935 601 signal = (x2>0.) ? par[k]*x2*exp(-x2+1.-z2) :0.0; // RCCR2
602 // signal =(x2>0.) ? par[k]*x2*x2*exp(-2*x2+2.-z2 ):0.0;//RCCR
603 }else if( electronics == 2 ) { // OLA
bf3f2830 604 x2 = (x-par[k+1])*(x-par[k+1])/t2;
50d05d7b 605 signal = par[k] * exp( -x2 - z2 );
606 } else {
aacedc3e 607 Warning("PeakFunc","Wrong SDD Electronics = %d",
608 electronics);
50d05d7b 609 // exit( 1 );
610 } // end if electronicx
24a1c341 611 spe[x*zdim+z] += signal;
612 if( integral != 0 ) integral[i] += signal;
42da2935 613 } // end for x
614 } // end for z
24a1c341 615 k += knParam;
42da2935 616 } // end for i
24a1c341 617 return;
a1f090e0 618}
42da2935 619//__________________________________________________________________________
aacedc3e 620Double_t AliITSClusterFinderSDD::ChiSqr( Int_t xdim, Int_t zdim, Double_t *spe,
621 Double_t *speFit ) const{
42da2935 622 // EVALUATES UNNORMALIZED CHI-SQUARED
aacedc3e 623 Double_t chi2 = 0.;
42da2935 624 for( Int_t z=0; z<zdim; z++ ){
50d05d7b 625 for( Int_t x=1; x<xdim-1; x++ ){
626 Int_t index = x*zdim+z;
aacedc3e 627 Double_t tmp = spe[index] - speFit[index];
50d05d7b 628 chi2 += tmp*tmp;
629 } // end for x
42da2935 630 } // end for z
631 return( chi2 );
a1f090e0 632}
42da2935 633//_______________________________________________________________________
aacedc3e 634void AliITSClusterFinderSDD::Minim( Int_t xdim, Int_t zdim, Double_t *param,
635 Double_t *prm0,Double_t *steprm,
636 Double_t *chisqr,Double_t *spe,
637 Double_t *speFit ){
42da2935 638 //
639 Int_t k, nnn, mmm, i;
aacedc3e 640 Double_t p1, delta, d1, chisq1, p2, chisq2, t, p3, chisq3, a, b, p0, chisqt;
42da2935 641 const Int_t knParam = 5;
642 Int_t npeak = (Int_t)param[0];
643 for( k=1; k<(npeak*knParam+1); k++ ) prm0[k] = param[k];
644 for( k=1; k<(npeak*knParam+1); k++ ){
50d05d7b 645 p1 = param[k];
646 delta = steprm[k];
647 d1 = delta;
648 // ENSURE THAT STEP SIZE IS SENSIBLY LARGER THAN MACHINE ROUND OFF
649 if( fabs( p1 ) > 1.0E-6 )
650 if ( fabs( delta/p1 ) < 1.0E-4 ) delta = p1/1000;
aacedc3e 651 else delta = (Double_t)1.0E-4;
50d05d7b 652 // EVALUATE CHI-SQUARED AT FIRST TWO SEARCH POINTS
653 PeakFunc( xdim, zdim, param, speFit );
654 chisq1 = ChiSqr( xdim, zdim, spe, speFit );
655 p2 = p1+delta;
656 param[k] = p2;
657 PeakFunc( xdim, zdim, param, speFit );
658 chisq2 = ChiSqr( xdim, zdim, spe, speFit );
659 if( chisq1 < chisq2 ){
660 // REVERSE DIRECTION OF SEARCH IF CHI-SQUARED IS INCREASING
661 delta = -delta;
662 t = p1;
663 p1 = p2;
664 p2 = t;
665 t = chisq1;
666 chisq1 = chisq2;
667 chisq2 = t;
668 } // end if
669 i = 1; nnn = 0;
670 do { // INCREMENT param(K) UNTIL CHI-SQUARED STARTS TO INCREASE
671 nnn++;
672 p3 = p2 + delta;
673 mmm = nnn - (nnn/5)*5; // multiplo de 5
674 if( mmm == 0 ){
675 d1 = delta;
676 // INCREASE STEP SIZE IF STEPPING TOWARDS MINIMUM IS TOO SLOW
677 delta *= 5;
678 } // end if
679 param[k] = p3;
680 // Constrain paramiters
681 Int_t kpos = (k-1) % knParam;
682 switch( kpos ){
683 case 0 :
684 if( param[k] <= 20 ) param[k] = fMinPeak;
685 break;
686 case 1 :
687 if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k];
688 break;
689 case 2 :
690 if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k];
691 break;
692 case 3 :
693 if( param[k] < .5 ) param[k] = .5;
694 break;
695 case 4 :
aacedc3e 696 if( param[k] < .288 ) param[k] = .288;// 1/sqrt(12) = 0.288
50d05d7b 697 if( param[k] > zdim*.5 ) param[k] = zdim*.5;
698 break;
699 }; // end switch
700 PeakFunc( xdim, zdim, param, speFit );
701 chisq3 = ChiSqr( xdim, zdim, spe, speFit );
702 if( chisq3 < chisq2 && nnn < 50 ){
703 p1 = p2;
704 p2 = p3;
705 chisq1 = chisq2;
706 chisq2 = chisq3;
707 }else i=0;
708 } while( i );
709 // FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS
710 a = chisq1*(p2-p3)+chisq2*(p3-p1)+chisq3*(p1-p2);
711 b = chisq1*(p2*p2-p3*p3)+chisq2*(p3*p3-p1*p1)+chisq3*(p1*p1-p2*p2);
aacedc3e 712 if( a!=0 ) p0 = (Double_t)(0.5*b/a);
50d05d7b 713 else p0 = 10000;
714 //--IN CASE OF NEARLY EQUAL CHI-SQUARED AND TOO SMALL STEP SIZE PREVENT
715 // ERRONEOUS EVALUATION OF PARABOLA MINIMUM
716 //---NEXT TWO LINES CAN BE OMITTED FOR HIGHER PRECISION MACHINES
aacedc3e 717 //dp = (Double_t) max (fabs(p3-p2), fabs(p2-p1));
50d05d7b 718 //if( fabs( p2-p0 ) > dp ) p0 = p2;
719 param[k] = p0;
720 // Constrain paramiters
721 Int_t kpos = (k-1) % knParam;
722 switch( kpos ){
723 case 0 :
724 if( param[k] <= 20 ) param[k] = fMinPeak;
725 break;
726 case 1 :
727 if( fabs( param[k] - prm0[k] ) > 1.5 ) param[k] = prm0[k];
728 break;
729 case 2 :
730 if( fabs( param[k] - prm0[k] ) > 1. ) param[k] = prm0[k];
731 break;
732 case 3 :
733 if( param[k] < .5 ) param[k] = .5;
734 break;
735 case 4 :
736 if( param[k] < .288 ) param[k] = .288; // 1/sqrt(12) = 0.288
737 if( param[k] > zdim*.5 ) param[k] = zdim*.5;
738 break;
739 }; // end switch
740 PeakFunc( xdim, zdim, param, speFit );
741 chisqt = ChiSqr( xdim, zdim, spe, speFit );
742 // DO NOT ALLOW ERRONEOUS INTERPOLATION
743 if( chisqt <= *chisqr ) *chisqr = chisqt;
744 else param[k] = prm0[k];
745 // OPTIMIZE SEARCH STEP FOR EVENTUAL NEXT CALL OF MINIM
746 steprm[k] = (param[k]-prm0[k])/5;
747 if( steprm[k] >= d1 ) steprm[k] = d1/5;
42da2935 748 } // end for k
749 // EVALUATE FIT AND CHI-SQUARED FOR OPTIMIZED PARAMETERS
750 PeakFunc( xdim, zdim, param, speFit );
751 *chisqr = ChiSqr( xdim, zdim, spe, speFit );
752 return;
a1f090e0 753}
42da2935 754//_________________________________________________________________________
755Int_t AliITSClusterFinderSDD::NoLinearFit( Int_t xdim, Int_t zdim,
aacedc3e 756 Double_t *param, Double_t *spe,
757 Int_t *niter, Double_t *chir ){
42da2935 758 // fit method from Comput. Phys. Commun 46(1987) 149
aacedc3e 759 const Double_t kchilmt = 0.01; // relative accuracy
50d05d7b 760 const Int_t knel = 3; // for parabolic minimization
761 const Int_t knstop = 50; // Max. iteration number
42da2935 762 const Int_t knParam = 5;
763 Int_t npeak = (Int_t)param[0];
764 // RETURN IF NUMBER OF DEGREES OF FREEDOM IS NOT POSITIVE
765 if( (xdim*zdim - npeak*knParam) <= 0 ) return( -1 );
aacedc3e 766 Double_t degFree = (xdim*zdim - npeak*knParam)-1;
42da2935 767 Int_t n, k, iterNum = 0;
aacedc3e 768 Double_t *prm0 = new Double_t[npeak*knParam+1];
769 Double_t *step = new Double_t[npeak*knParam+1];
770 Double_t *schi = new Double_t[npeak*knParam+1];
771 Double_t *sprm[3];
772 sprm[0] = new Double_t[npeak*knParam+1];
773 sprm[1] = new Double_t[npeak*knParam+1];
774 sprm[2] = new Double_t[npeak*knParam+1];
775 Double_t chi0, chi1, reldif, a, b, prmin, dp;
776 Double_t *speFit = new Double_t[ xdim*zdim ];
42da2935 777 PeakFunc( xdim, zdim, param, speFit );
778 chi0 = ChiSqr( xdim, zdim, spe, speFit );
779 chi1 = chi0;
780 for( k=1; k<(npeak*knParam+1); k++) prm0[k] = param[k];
50d05d7b 781 for( k=1 ; k<(npeak*knParam+1); k+=knParam ){
782 step[k] = param[k] / 20.0 ;
783 step[k+1] = param[k+1] / 50.0;
784 step[k+2] = param[k+2] / 50.0;
785 step[k+3] = param[k+3] / 20.0;
786 step[k+4] = param[k+4] / 20.0;
787 } // end for k
42da2935 788 Int_t out = 0;
789 do{
50d05d7b 790 iterNum++;
791 chi0 = chi1;
792 Minim( xdim, zdim, param, prm0, step, &chi1, spe, speFit );
aacedc3e 793 reldif = ( chi1 > 0 ) ? ((Double_t) fabs( chi1-chi0)/chi1 ) : 0;
50d05d7b 794 // EXIT conditions
795 if( reldif < (float) kchilmt ){
796 *chir = (chi1>0) ? (float) TMath::Sqrt (chi1/degFree) :0;
797 *niter = iterNum;
798 out = 0;
799 break;
800 } // end if
801 if( (reldif < (float)(5*kchilmt)) && (iterNum > knstop) ){
802 *chir = (chi1>0) ?(float) TMath::Sqrt (chi1/degFree):0;
803 *niter = iterNum;
804 out = 0;
805 break;
806 } // end if
807 if( iterNum > 5*knstop ){
808 *chir = (chi1>0) ?(float) TMath::Sqrt (chi1/degFree):0;
809 *niter = iterNum;
810 out = 1;
811 break;
812 } // end if
813 if( iterNum <= knel ) continue;
814 n = iterNum - (iterNum/knel)*knel; // EXTRAPOLATION LIMIT COUNTER N
815 if( n > 3 || n == 0 ) continue;
816 schi[n-1] = chi1;
817 for( k=1; k<(npeak*knParam+1); k++ ) sprm[n-1][k] = param[k];
818 if( n != 3 ) continue;
819 // -EVALUATE EXTRAPOLATED VALUE OF EACH PARAMETER BY FINDING MINIMUM OF
820 // PARABOLA DEFINED BY LAST THREE CALLS OF MINIM
821 for( k=1; k<(npeak*knParam+1); k++ ){
aacedc3e 822 Double_t tmp0 = sprm[0][k];
823 Double_t tmp1 = sprm[1][k];
824 Double_t tmp2 = sprm[2][k];
50d05d7b 825 a = schi[0]*(tmp1-tmp2) + schi[1]*(tmp2-tmp0);
826 a += (schi[2]*(tmp0-tmp1));
827 b = schi[0]*(tmp1*tmp1-tmp2*tmp2);
828 b += (schi[1]*(tmp2*tmp2-tmp0*tmp0)+(schi[2]*
829 (tmp0*tmp0-tmp1*tmp1)));
830 if ((double)a < 1.0E-6) prmin = 0;
831 else prmin = (float) (0.5*b/a);
832 dp = 5*(tmp2-tmp0);
833 if( fabs(prmin-tmp2) > fabs(dp) ) prmin = tmp2+dp;
834 param[k] = prmin;
835 step[k] = dp/10; // OPTIMIZE SEARCH STEP
836 } // end for k
42da2935 837 } while( kTRUE );
838 delete [] prm0;
839 delete [] step;
840 delete [] schi;
841 delete [] sprm[0];
842 delete [] sprm[1];
843 delete [] sprm[2];
844 delete [] speFit;
845 return( out );
a1f090e0 846}
50d05d7b 847
42da2935 848//______________________________________________________________________
04366a57 849void AliITSClusterFinderSDD::ResolveClusters(){
42da2935 850 // The function to resolve clusters if the clusters overlapping exists
24a1c341 851 Int_t i;
42da2935 852 // get number of clusters for this module
aacedc3e 853 Int_t nofClusters = NClusters();
42da2935 854 nofClusters -= fNclusters;
aacedc3e 855 Int_t fNofMaps = GetSeg()->Npz();
42da2935 856 Int_t fNofAnodes = fNofMaps/2;
aacedc3e 857 //Int_t fMaxNofSamples = GetSeg()->Npx();
42da2935 858 Int_t dummy=0;
aacedc3e 859 Double_t fTimeStep = GetSeg()->Dpx( dummy );
860 Double_t fSddLength = GetSeg()->Dx();
861 Double_t fDriftSpeed = GetResp()->DriftSpeed();
862 Double_t anodePitch = GetSeg()->Dpz( dummy );
863 Double_t n, baseline;
864 GetResp()->GetNoiseParam( n, baseline );
865 Int_t electronics = GetResp()->Electronics(); // 1 = PASCAL, 2 = OLA
50d05d7b 866
42da2935 867 for( Int_t j=0; j<nofClusters; j++ ){
50d05d7b 868 // get cluster information
aacedc3e 869 AliITSRawClusterSDD *clusterJ=(AliITSRawClusterSDD*) Cluster(j);
50d05d7b 870 Int_t astart = clusterJ->Astart();
871 Int_t astop = clusterJ->Astop();
872 Int_t tstart = clusterJ->Tstartf();
873 Int_t tstop = clusterJ->Tstopf();
874 Int_t wing = (Int_t)clusterJ->W();
875 if( wing == 2 ){
876 astart += fNofAnodes;
877 astop += fNofAnodes;
878 } // end if
879 Int_t xdim = tstop-tstart+3;
880 Int_t zdim = astop-astart+3;
d86f531c 881 if( xdim > 50 || zdim > 30 ) {
04366a57 882 Warning("ResolveClusters","xdim: %d , zdim: %d ",xdim,zdim);
d86f531c 883 continue;
884 }
aacedc3e 885 Double_t *sp = new Double_t[ xdim*zdim+1 ];
886 memset( sp, 0, sizeof(Double_t)*(xdim*zdim+1) );
50d05d7b 887
888 // make a local map from cluster region
889 for( Int_t ianode=astart; ianode<=astop; ianode++ ){
890 for( Int_t itime=tstart; itime<=tstop; itime++ ){
aacedc3e 891 Double_t fadc = Map()->GetSignal( ianode, itime );
50d05d7b 892 if( fadc > baseline ) fadc -= (Double_t)baseline;
893 else fadc = 0.;
894 Int_t index = (itime-tstart+1)*zdim+(ianode-astart+1);
895 sp[index] = fadc;
896 } // time loop
897 } // anode loop
898
899 // search peaks on cluster
900 const Int_t kNp = 150;
901 Int_t peakX1[kNp];
902 Int_t peakZ1[kNp];
aacedc3e 903 Double_t peakAmp1[kNp];
50d05d7b 904 Int_t npeak = SearchPeak(sp,xdim,zdim,peakX1,peakZ1,peakAmp1,fMinPeak);
905
906 // if multiple peaks, split cluster
aacedc3e 907 if( npeak >= 1 ){
50d05d7b 908 // cout << "npeak " << npeak << endl;
909 // clusterJ->PrintInfo();
aacedc3e 910 Double_t *par = new Double_t[npeak*5+1];
911 par[0] = (Double_t)npeak;
48058160 912 // Initial parameters in cell dimentions
50d05d7b 913 Int_t k1 = 1;
914 for( i=0; i<npeak; i++ ){
915 par[k1] = peakAmp1[i];
916 par[k1+1] = peakX1[i]; // local time pos. [timebin]
917 par[k1+2] = peakZ1[i]; // local anode pos. [anodepitch]
aacedc3e 918 if( electronics == 1 ) par[k1+3] = 2.; // PASCAL
919 else if(electronics==2) par[k1+3] = 0.7;//tau [timebin] OLA
50d05d7b 920 par[k1+4] = .4; // sigma [anodepich]
aacedc3e 921 k1 += 5;
50d05d7b 922 } // end for i
923 Int_t niter;
aacedc3e 924 Double_t chir;
50d05d7b 925 NoLinearFit( xdim, zdim, par, sp, &niter, &chir );
aacedc3e 926 Double_t peakX[kNp];
927 Double_t peakZ[kNp];
928 Double_t sigma[kNp];
929 Double_t tau[kNp];
930 Double_t peakAmp[kNp];
931 Double_t integral[kNp];
50d05d7b 932 //get integrals => charge for each peak
933 PeakFunc( xdim, zdim, par, sp, integral );
934 k1 = 1;
935 for( i=0; i<npeak; i++ ){
936 peakAmp[i] = par[k1];
d86f531c 937 peakX[i] = par[k1+1];
938 peakZ[i] = par[k1+2];
939 tau[i] = par[k1+3];
940 sigma[i] = par[k1+4];
50d05d7b 941 k1+=5;
942 } // end for i
943 // calculate parameter for new clusters
944 for( i=0; i<npeak; i++ ){
945 AliITSRawClusterSDD clusterI( *clusterJ );
d86f531c 946
50d05d7b 947 Int_t newAnode = peakZ1[i]-1 + astart;
d86f531c 948
949 // Int_t newiTime = peakX1[i]-1 + tstart;
950 // Int_t shift = (Int_t)(fTimeCorr/fTimeStep + 0.5);
951 // if( newiTime > shift && newiTime < (fMaxNofSamples-shift) )
952 // shift = 0;
aacedc3e 953 // Int_t peakpos = Map()->GetHitIndex(newAnode,newiTime+shift );
d86f531c 954 // clusterI.SetPeakPos( peakpos );
955
50d05d7b 956 clusterI.SetPeakAmpl( peakAmp1[i] );
aacedc3e 957 Double_t newAnodef = peakZ[i] - 0.5 + astart;
958 Double_t newiTimef = peakX[i] - 1 + tstart;
50d05d7b 959 if( wing == 2 ) newAnodef -= fNofAnodes;
aacedc3e 960 Double_t anodePath = (newAnodef - fNofAnodes/2)*anodePitch;
50d05d7b 961 newiTimef *= fTimeStep;
962 if( newiTimef > fTimeCorr ) newiTimef -= fTimeCorr;
963 if( electronics == 1 ){
48058160 964 // newiTimef *= 0.999438; // PASCAL
965 // newiTimef += (6./fDriftSpeed - newiTimef/3000.);
50d05d7b 966 }else if( electronics == 2 )
967 newiTimef *= 0.99714; // OLA
d86f531c 968
969 Int_t timeBin = (Int_t)(newiTimef/fTimeStep+0.5);
aacedc3e 970 Int_t peakpos = Map()->GetHitIndex( newAnode, timeBin );
d86f531c 971 if( peakpos < 0 ) {
972 for( Int_t ii=0; ii<3; ii++ ) {
aacedc3e 973 peakpos = Map()->GetHitIndex( newAnode, timeBin+ii );
d86f531c 974 if( peakpos > 0 ) break;
aacedc3e 975 peakpos = Map()->GetHitIndex( newAnode, timeBin-ii );
d86f531c 976 if( peakpos > 0 ) break;
977 }
978 }
979
980 if( peakpos < 0 ) {
04366a57 981 //Warning("ResolveClusters",
aacedc3e 982 // "Digit not found for cluster");
983 //if(GetDebug(3)) clusterI.PrintInfo();
d86f531c 984 continue;
985 }
986 clusterI.SetPeakPos( peakpos );
aacedc3e 987 Double_t driftPath = fSddLength - newiTimef * fDriftSpeed;
988 Double_t sign = ( wing == 1 ) ? -1. : 1.;
50d05d7b 989 clusterI.SetX( driftPath*sign * 0.0001 );
990 clusterI.SetZ( anodePath * 0.0001 );
991 clusterI.SetAnode( newAnodef );
992 clusterI.SetTime( newiTimef );
993 clusterI.SetAsigma( sigma[i]*anodePitch );
994 clusterI.SetTsigma( tau[i]*fTimeStep );
995 clusterI.SetQ( integral[i] );
d86f531c 996
04366a57 997 fITS->AddCluster( 1, &clusterI );
50d05d7b 998 } // end for i
aacedc3e 999 Clusters()->RemoveAt( j );
50d05d7b 1000 delete [] par;
48058160 1001 } else { // something odd
04366a57 1002 Warning( "ResolveClusters",
aacedc3e 1003 "--- Peak not found!!!! minpeak=%d ,cluster peak= %f"
1004 " , module= %d",
1005 fMinPeak, clusterJ->PeakAmpl(),GetModule());
d86f531c 1006 clusterJ->PrintInfo();
04366a57 1007 Warning( "ResolveClusters"," xdim= %d zdim= %d", xdim-2, zdim-2 );
50d05d7b 1008 }
1009 delete [] sp;
42da2935 1010 } // cluster loop
aacedc3e 1011 Clusters()->Compress();
1012// Map()->ClearMap();
a1f090e0 1013}
42da2935 1014//________________________________________________________________________
1015void AliITSClusterFinderSDD::GroupClusters(){
1016 // group clusters
1017 Int_t dummy=0;
aacedc3e 1018 Double_t fTimeStep = GetSeg()->Dpx(dummy);
42da2935 1019 // get number of clusters for this module
aacedc3e 1020 Int_t nofClusters = NClusters();
42da2935 1021 nofClusters -= fNclusters;
1022 AliITSRawClusterSDD *clusterI;
1023 AliITSRawClusterSDD *clusterJ;
1024 Int_t *label = new Int_t [nofClusters];
1025 Int_t i,j;
1026 for(i=0; i<nofClusters; i++) label[i] = 0;
1027 for(i=0; i<nofClusters; i++) {
50d05d7b 1028 if(label[i] != 0) continue;
1029 for(j=i+1; j<nofClusters; j++) {
1030 if(label[j] != 0) continue;
aacedc3e 1031 clusterI = (AliITSRawClusterSDD*) Cluster(i);
1032 clusterJ = (AliITSRawClusterSDD*) Cluster(j);
50d05d7b 1033 // 1.3 good
1034 if(clusterI->T() < fTimeStep*60) fDAnode = 4.2; // TB 3.2
1035 if(clusterI->T() < fTimeStep*10) fDAnode = 1.5; // TB 1.
1036 Bool_t pair = clusterI->Brother(clusterJ,fDAnode,fDTime);
1037 if(!pair) continue;
aacedc3e 1038 if(GetDebug(4)){
1039 clusterI->PrintInfo();
1040 clusterJ->PrintInfo();
1041 } // end if GetDebug
50d05d7b 1042 clusterI->Add(clusterJ);
1043 label[j] = 1;
aacedc3e 1044 Clusters()->RemoveAt(j);
50d05d7b 1045 j=i; // <- Ernesto
1046 } // J clusters
1047 label[i] = 1;
42da2935 1048 } // I clusters
aacedc3e 1049 Clusters()->Compress();
42da2935 1050
1051 delete [] label;
1052 return;
b0f5e3fc 1053}
42da2935 1054//________________________________________________________________________
1055void AliITSClusterFinderSDD::SelectClusters(){
1056 // get number of clusters for this module
aacedc3e 1057 Int_t nofClusters = NClusters();
b0f5e3fc 1058
42da2935 1059 nofClusters -= fNclusters;
1060 Int_t i;
1061 for(i=0; i<nofClusters; i++) {
aacedc3e 1062 AliITSRawClusterSDD *clusterI =(AliITSRawClusterSDD*) Cluster(i);
50d05d7b 1063 Int_t rmflg = 0;
aacedc3e 1064 Double_t wy = 0.;
50d05d7b 1065 if(clusterI->Anodes() != 0.) {
aacedc3e 1066 wy = ((Double_t) clusterI->Samples())/clusterI->Anodes();
50d05d7b 1067 } // end if
1068 Int_t amp = (Int_t) clusterI->PeakAmpl();
1069 Int_t cha = (Int_t) clusterI->Q();
1070 if(amp < fMinPeak) rmflg = 1;
1071 if(cha < fMinCharge) rmflg = 1;
1072 if(wy < fMinNCells) rmflg = 1;
1073 //if(wy > fMaxNCells) rmflg = 1;
aacedc3e 1074 if(rmflg) Clusters()->RemoveAt(i);
42da2935 1075 } // I clusters
aacedc3e 1076 Clusters()->Compress();
42da2935 1077 return;
b0f5e3fc 1078}
42da2935 1079
42da2935 1080//______________________________________________________________________
1081void AliITSClusterFinderSDD::GetRecPoints(){
1082 // get rec points
04366a57 1083
42da2935 1084 // get number of clusters for this module
aacedc3e 1085 Int_t nofClusters = NClusters();
42da2935 1086 nofClusters -= fNclusters;
aacedc3e 1087 const Double_t kconvGeV = 1.e-6; // GeV -> KeV
1088 const Double_t kconv = 1.0e-4;
1089 const Double_t kRMSx = 38.0*kconv; // microns->cm ITS TDR Table 1.3
1090 const Double_t kRMSz = 28.0*kconv; // microns->cm ITS TDR Table 1.3
d86f531c 1091 Int_t i;
42da2935 1092 Int_t ix, iz, idx=-1;
1093 AliITSdigitSDD *dig=0;
aacedc3e 1094 Int_t ndigits=NDigits();
42da2935 1095 for(i=0; i<nofClusters; i++) {
aacedc3e 1096 AliITSRawClusterSDD *clusterI = (AliITSRawClusterSDD*)Cluster(i);
50d05d7b 1097 if(!clusterI) Error("SDD: GetRecPoints","i clusterI ",i,clusterI);
1098 if(clusterI) idx=clusterI->PeakPos();
1099 if(idx>ndigits) Error("SDD: GetRecPoints","idx ndigits",idx,ndigits);
1100 // try peak neighbours - to be done
aacedc3e 1101 if(idx&&idx<= ndigits) dig =(AliITSdigitSDD*)GetDigit(idx);
d86f531c 1102 if(!dig) {
50d05d7b 1103 // try cog
aacedc3e 1104 GetSeg()->GetPadIxz(clusterI->X(),clusterI->Z(),ix,iz);
1105 dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix-1);
50d05d7b 1106 // if null try neighbours
aacedc3e 1107 if (!dig) dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix);
1108 if (!dig) dig = (AliITSdigitSDD*)Map()->GetHit(iz-1,ix+1);
50d05d7b 1109 if (!dig) printf("SDD: cannot assign the track number!\n");
1110 } // end if !dig
1111 AliITSRecPoint rnew;
1112 rnew.SetX(clusterI->X());
1113 rnew.SetZ(clusterI->Z());
1114 rnew.SetQ(clusterI->Q()); // in KeV - should be ADC
1115 rnew.SetdEdX(kconvGeV*clusterI->Q());
1116 rnew.SetSigmaX2(kRMSx*kRMSx);
1117 rnew.SetSigmaZ2(kRMSz*kRMSz);
d86f531c 1118
ecee53fc 1119 if(dig) rnew.fTracks[0]=dig->GetTrack(0);
1120 if(dig) rnew.fTracks[1]=dig->GetTrack(1);
1121 if(dig) rnew.fTracks[2]=dig->GetTrack(2);
507cf1a7 1122
04366a57 1123 fITS->AddRecPoint(rnew);
42da2935 1124 } // I clusters
aacedc3e 1125// Map()->ClearMap();
b0f5e3fc 1126}
42da2935 1127//______________________________________________________________________
1128void AliITSClusterFinderSDD::FindRawClusters(Int_t mod){
1129 // find raw clusters
50d05d7b 1130
aacedc3e 1131 SetModule(mod);
a1f090e0 1132 Find1DClustersE();
b0f5e3fc 1133 GroupClusters();
1134 SelectClusters();
04366a57 1135 ResolveClusters();
b0f5e3fc 1136 GetRecPoints();
1137}
42da2935 1138//_______________________________________________________________________
bf3f2830 1139void AliITSClusterFinderSDD::Print() const{
42da2935 1140 // Print SDD cluster finder Parameters
1141
1142 cout << "**************************************************" << endl;
1143 cout << " Silicon Drift Detector Cluster Finder Parameters " << endl;
1144 cout << "**************************************************" << endl;
1145 cout << "Number of Clusters: " << fNclusters << endl;
1146 cout << "Anode Tolerance: " << fDAnode << endl;
1147 cout << "Time Tolerance: " << fDTime << endl;
1148 cout << "Time correction (electronics): " << fTimeCorr << endl;
1149 cout << "Cut Amplitude (threshold): " << fCutAmplitude << endl;
1150 cout << "Minimum Amplitude: " << fMinPeak << endl;
1151 cout << "Minimum Charge: " << fMinCharge << endl;
1152 cout << "Minimum number of cells/clusters: " << fMinNCells << endl;
1153 cout << "Maximum number of cells/clusters: " << fMaxNCells << endl;
1154 cout << "**************************************************" << endl;
a1f090e0 1155}