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