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