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