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