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