Changed to Copy(TObject &A) to appease the HP-compiler
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
CommitLineData
f7336fa3 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 **************************************************************************/
15
16/*
17$Log$
43da34c0 18Revision 1.6 2000/06/09 11:10:07 cblume
19Compiler warnings and coding conventions, next round
20
dd9a6ee3 21Revision 1.5 2000/06/08 18:32:58 cblume
22Make code compliant to coding conventions
23
8230f242 24Revision 1.4 2000/06/07 16:27:01 cblume
25Try to remove compiler warnings on Sun and HP
26
9d0b222b 27Revision 1.3 2000/05/08 16:17:27 cblume
28Merge TRD-develop
29
6f1e466d 30Revision 1.1.4.1 2000/05/08 15:09:01 cblume
31Introduce AliTRDdigitsManager
32
c0dd96c3 33Revision 1.1 2000/02/28 18:58:54 cblume
34Add new TRD classes
35
f7336fa3 36*/
37
38///////////////////////////////////////////////////////////////////////////////
39// //
40// TRD cluster finder for the slow simulator.
41// //
42///////////////////////////////////////////////////////////////////////////////
43
44#include <TF1.h>
45
46#include "AliTRDclusterizerV1.h"
47#include "AliTRDmatrix.h"
48#include "AliTRDgeometry.h"
49#include "AliTRDdigitizer.h"
50#include "AliTRDrecPoint.h"
6f1e466d 51#include "AliTRDdataArrayF.h"
f7336fa3 52
53ClassImp(AliTRDclusterizerV1)
54
55//_____________________________________________________________________________
56AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
57{
58 //
59 // AliTRDclusterizerV1 default constructor
60 //
61
6f1e466d 62 fDigitsManager = NULL;
f7336fa3 63
64}
65
66//_____________________________________________________________________________
67AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
68 :AliTRDclusterizer(name,title)
69{
70 //
71 // AliTRDclusterizerV1 default constructor
72 //
73
6f1e466d 74 fDigitsManager = new AliTRDdigitsManager();
f7336fa3 75
76 Init();
77
78}
79
8230f242 80//_____________________________________________________________________________
dd9a6ee3 81AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
8230f242 82{
83 //
84 // AliTRDclusterizerV1 copy constructor
85 //
86
dd9a6ee3 87 ((AliTRDclusterizerV1 &) c).Copy(*this);
8230f242 88
89}
90
f7336fa3 91//_____________________________________________________________________________
92AliTRDclusterizerV1::~AliTRDclusterizerV1()
93{
8230f242 94 //
95 // AliTRDclusterizerV1 destructor
96 //
f7336fa3 97
6f1e466d 98 if (fDigitsManager) {
99 delete fDigitsManager;
f7336fa3 100 }
101
102}
103
dd9a6ee3 104//_____________________________________________________________________________
105AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
106{
107 //
108 // Assignment operator
109 //
110
111 if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
112 return *this;
113
114}
115
8230f242 116//_____________________________________________________________________________
43da34c0 117void AliTRDclusterizerV1::Copy(TObject &c)
8230f242 118{
119 //
120 // Copy function
121 //
122
43da34c0 123 ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
124 ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
125 ((AliTRDclusterizerV1 &) c).fClusMethod = fClusMethod;
126 ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
8230f242 127
128 AliTRDclusterizer::Copy(c);
129
130}
131
f7336fa3 132//_____________________________________________________________________________
133void AliTRDclusterizerV1::Init()
134{
135 //
136 // Initializes the cluster finder
137 //
138
139 // The default parameter for the clustering
140 fClusMaxThresh = 5.0;
141 fClusSigThresh = 2.0;
142 fClusMethod = 1;
143
144}
145
146//_____________________________________________________________________________
147Bool_t AliTRDclusterizerV1::ReadDigits()
148{
149 //
150 // Reads the digits arrays from the input aliroot file
151 //
152
153 if (!fInputFile) {
154 printf("AliTRDclusterizerV1::ReadDigits -- ");
155 printf("No input file open\n");
156 return kFALSE;
157 }
158
f7336fa3 159 // Read in the digit arrays
6f1e466d 160 return (fDigitsManager->ReadDigits());
f7336fa3 161
162}
163
164//_____________________________________________________________________________
165Bool_t AliTRDclusterizerV1::MakeCluster()
166{
167 //
168 // Generates the cluster.
169 //
170
171 Int_t row, col, time;
172
173 // Get the pointer to the detector class and check for version 1
8230f242 174 AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD");
175 if (trd->IsVersion() != 1) {
f7336fa3 176 printf("AliTRDclusterizerV1::MakeCluster -- ");
177 printf("TRD must be version 1 (slow simulator).\n");
178 return kFALSE;
179 }
180
181 // Get the geometry
8230f242 182 AliTRDgeometry *geo = trd->GetGeometry();
f7336fa3 183
184 printf("AliTRDclusterizerV1::MakeCluster -- ");
185 printf("Start creating clusters.\n");
186
8230f242 187 AliTRDdataArrayI *digits;
f7336fa3 188
189 // Parameters
190 Float_t maxThresh = fClusMaxThresh; // threshold value for maximum
191 Float_t signalThresh = fClusSigThresh; // threshold value for digit signal
192 Int_t clusteringMethod = fClusMethod; // clustering method option (for testing)
193
194 // Iteration limit for unfolding procedure
8230f242 195 const Float_t kEpsilon = 0.01;
f7336fa3 196
8230f242 197 const Int_t kNclus = 3;
198 const Int_t kNsig = 5;
f7336fa3 199
200 Int_t chamBeg = 0;
201 Int_t chamEnd = kNcham;
8230f242 202 if (trd->GetSensChamber() >= 0) {
203 chamBeg = trd->GetSensChamber();
6f1e466d 204 chamEnd = chamBeg + 1;
f7336fa3 205 }
206 Int_t planBeg = 0;
207 Int_t planEnd = kNplan;
8230f242 208 if (trd->GetSensPlane() >= 0) {
209 planBeg = trd->GetSensPlane();
f7336fa3 210 planEnd = planBeg + 1;
211 }
212 Int_t sectBeg = 0;
213 Int_t sectEnd = kNsect;
f7336fa3 214
215 // *** Start clustering *** in every chamber
216 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
217 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
218 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
219
8230f242 220 if (trd->GetSensSector() >= 0) {
221 Int_t sens1 = trd->GetSensSector();
222 Int_t sens2 = sens1 + trd->GetSensSectorRange();
9d0b222b 223 sens2 -= ((Int_t) (sens2 / kNsect)) * kNsect;
dd9a6ee3 224 if (sens1 < sens2) {
9d0b222b 225 if ((isect < sens1) || (isect >= sens2)) continue;
dd9a6ee3 226 }
227 else {
9d0b222b 228 if ((isect < sens1) && (isect >= sens2)) continue;
dd9a6ee3 229 }
9d0b222b 230 }
231
8230f242 232 Int_t idet = geo->GetDetector(iplan,icham,isect);
f7336fa3 233
234 Int_t nClusters = 0;
235 printf("AliTRDclusterizerV1::MakeCluster -- ");
236 printf("Analyzing chamber %d, plane %d, sector %d.\n"
237 ,icham,iplan,isect);
238
8230f242 239 Int_t nRowMax = geo->GetRowMax(iplan,icham,isect);
240 Int_t nColMax = geo->GetColMax(iplan);
241 Int_t nTimeMax = geo->GetTimeMax();
f7336fa3 242
243 // Create a detector matrix to keep maxima
244 AliTRDmatrix *digitMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
245 ,isect,icham,iplan);
246 // Create a matrix to contain maximum flags
247 AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax
248 ,isect,icham,iplan);
249
250 // Read in the digits
8230f242 251 digits = fDigitsManager->GetDigits(idet);
f7336fa3 252
253 // Loop through the detector pixel
254 for (time = 0; time < nTimeMax; time++) {
255 for ( col = 0; col < nColMax; col++) {
256 for ( row = 0; row < nRowMax; row++) {
257
8230f242 258 Int_t signal = digits->GetData(row,col,time);
259 Int_t index = digits->GetIndex(row,col,time);
f7336fa3 260
261 // Fill the detector matrix
262 if (signal > signalThresh) {
263 // Store the signal amplitude
264 digitMatrix->SetSignal(row,col,time,signal);
265 // Store the digits number
266 digitMatrix->AddTrack(row,col,time,index);
267 }
268
269 }
270 }
271 }
272
273 // Loop chamber and find maxima in digitMatrix
274 for ( row = 0; row < nRowMax; row++) {
275 for ( col = 1; col < nColMax; col++) {
276 for (time = 0; time < nTimeMax; time++) {
277
278 if (digitMatrix->GetSignal(row,col,time)
279 < digitMatrix->GetSignal(row,col - 1,time)) {
280 // really maximum?
281 if (col > 1) {
282 if (digitMatrix->GetSignal(row,col - 2,time)
283 < digitMatrix->GetSignal(row,col - 1,time)) {
284 // yes, so set maximum flag
285 maximaMatrix->SetSignal(row,col - 1,time,1);
286 }
287 else maximaMatrix->SetSignal(row,col - 1,time,0);
288 }
289 }
290
291 } // time
292 } // col
293 } // row
294
295 // now check maxima and calculate cluster position
296 for ( row = 0; row < nRowMax; row++) {
297 for ( col = 1; col < nColMax; col++) {
298 for (time = 0; time < nTimeMax; time++) {
299
300 if ((maximaMatrix->GetSignal(row,col,time) > 0)
301 && (digitMatrix->GetSignal(row,col,time) > maxThresh)) {
302
303 // Ratio resulting from unfolding
8230f242 304 Float_t ratio = 0;
f7336fa3 305 // Signals on max and neighbouring pads
8230f242 306 Float_t padSignal[kNsig] = {0};
f7336fa3 307 // Signals from cluster
8230f242 308 Float_t clusterSignal[kNclus] = {0};
f7336fa3 309 // Cluster pad info
8230f242 310 Float_t clusterPads[kNclus] = {0};
f7336fa3 311 // Cluster digit info
8230f242 312 Int_t clusterDigit[kNclus] = {0};
f7336fa3 313
9d0b222b 314 Int_t iPad;
8230f242 315 for (iPad = 0; iPad < kNclus; iPad++) {
f7336fa3 316 clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
317 clusterDigit[iPad] = digitMatrix->GetTrack(row,col-1+iPad,time,0);
318 }
319
320 // neighbouring maximum on right side?
321 if (col < nColMax - 2) {
322 if (maximaMatrix->GetSignal(row,col + 2,time) > 0) {
323
9d0b222b 324 for (iPad = 0; iPad < 5; iPad++) {
f7336fa3 325 padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
326 }
327
328 // unfold:
8230f242 329 ratio = Unfold(kEpsilon, padSignal);
f7336fa3 330
331 // set signal on overlapping pad to ratio
332 clusterSignal[2] *= ratio;
333
334 }
335 }
336
337 // Calculate the position of the cluster
338 switch (clusteringMethod) {
339 case 1:
340 // method 1: simply center of mass
341 clusterPads[0] = row + 0.5;
342 clusterPads[1] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) /
c0dd96c3 343 (clusterSignal[0] + clusterSignal[1] + clusterSignal[2]);
f7336fa3 344 clusterPads[2] = time + 0.5;
345
346 nClusters++;
347 break;
348 case 2:
349 // method 2: integral gauss fit on 3 pads
350 TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads"
351 , 5, -1.5, 3.5);
352 for (Int_t iCol = -1; iCol <= 3; iCol++) {
353 if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1;
354 hPadCharges->Fill(iCol, clusterSignal[iCol]);
355 }
356 hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5);
357 TF1 *fPadChargeFit = hPadCharges->GetFunction("gaus");
358 Double_t colMean = fPadChargeFit->GetParameter(1);
359
360 clusterPads[0] = row + 0.5;
361 clusterPads[1] = col - 1.5 + colMean;
362 clusterPads[2] = time + 0.5;
363
364 delete hPadCharges;
365
366 nClusters++;
367 break;
368 }
369
370 Float_t clusterCharge = clusterSignal[0]
371 + clusterSignal[1]
372 + clusterSignal[2];
373
374 // Add the cluster to the output array
8230f242 375 trd->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge);
f7336fa3 376
377 }
378 } // time
379 } // col
380 } // row
381
382 printf("AliTRDclusterizerV1::MakeCluster -- ");
383 printf("Number of clusters found: %d\n",nClusters);
384
385 delete digitMatrix;
386 delete maximaMatrix;
387
388 } // isect
389 } // iplan
390 } // icham
391
392 printf("AliTRDclusterizerV1::MakeCluster -- ");
393 printf("Total number of points found: %d\n"
8230f242 394 ,trd->RecPoints()->GetEntries());
f7336fa3 395
396 // Get the pointer to the cluster branch
8230f242 397 TTree *clusterTree = gAlice->TreeR();
f7336fa3 398
399 // Fill the cluster-branch
400 printf("AliTRDclusterizerV1::MakeCluster -- ");
401 printf("Fill the cluster tree.\n");
8230f242 402 clusterTree->Fill();
f7336fa3 403 printf("AliTRDclusterizerV1::MakeCluster -- ");
404 printf("Done.\n");
405
406 return kTRUE;
407
408}
409
410//_____________________________________________________________________________
411Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
412{
413 //
414 // Method to unfold neighbouring maxima.
415 // The charge ratio on the overlapping pad is calculated
416 // until there is no more change within the range given by eps.
417 // The resulting ratio is then returned to the calling method.
418 //
419
420 Int_t itStep = 0; // count iteration steps
421
422 Float_t ratio = 0.5; // start value for ratio
423 Float_t prevRatio = 0; // store previous ratio
424
425 Float_t newLeftSignal[3] = {0}; // array to store left cluster signal
426 Float_t newRightSignal[3] = {0}; // array to store right cluster signal
427
428 // start iteration:
429 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
430
431 itStep++;
432 prevRatio = ratio;
433
434 // cluster position according to charge ratio
435 Float_t maxLeft = (ratio*padSignal[2] - padSignal[0]) /
436 (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
437 Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) /
438 ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
439
440 // set cluster charge ratio
441 Float_t ampLeft = padSignal[1];
442 Float_t ampRight = padSignal[3];
443
444 // apply pad response to parameters
445 newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft);
446 newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft);
447 newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft);
448
449 newRightSignal[0] = ampRight*PadResponse(-1 - maxRight);
450 newRightSignal[1] = ampRight*PadResponse( 0 - maxRight);
451 newRightSignal[2] = ampRight*PadResponse( 1 - maxRight);
452
453 // calculate new overlapping ratio
454 ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]);
455
456 }
457
458 return ratio;
459
460}
461
462//_____________________________________________________________________________
463Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
464{
465 //
466 // The pad response for the chevron pads.
467 // We use a simple Gaussian approximation which should be good
468 // enough for our purpose.
469 //
470
471 // The parameters for the response function
8230f242 472 const Float_t kA = 0.8872;
473 const Float_t kB = -0.00573;
474 const Float_t kC = 0.454;
475 const Float_t kC2 = kC*kC;
f7336fa3 476
8230f242 477 Float_t pr = kA * (kB + TMath::Exp(-x*x / (2. * kC2)));
f7336fa3 478
479 return (pr);
480
481}