Lost changes from revision 1.13 recovered.
[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$
3e1a3ad8 18Revision 1.9 2000/11/01 14:53:20 cblume
19Merge with TRD-develop
20
793ff80c 21Revision 1.1.4.5 2000/10/15 23:40:01 cblume
22Remove AliTRDconst
23
24Revision 1.1.4.4 2000/10/06 16:49:46 cblume
25Made Getters const
26
27Revision 1.1.4.3 2000/10/04 16:34:58 cblume
28Replace include files by forward declarations
29
30Revision 1.1.4.2 2000/09/22 14:49:49 cblume
31Adapted to tracking code
32
33Revision 1.8 2000/10/02 21:28:19 fca
34Removal of useless dependecies via forward declarations
35
94de3818 36Revision 1.7 2000/06/27 13:08:50 cblume
37Changed to Copy(TObject &A) to appease the HP-compiler
38
43da34c0 39Revision 1.6 2000/06/09 11:10:07 cblume
40Compiler warnings and coding conventions, next round
41
dd9a6ee3 42Revision 1.5 2000/06/08 18:32:58 cblume
43Make code compliant to coding conventions
44
8230f242 45Revision 1.4 2000/06/07 16:27:01 cblume
46Try to remove compiler warnings on Sun and HP
47
9d0b222b 48Revision 1.3 2000/05/08 16:17:27 cblume
49Merge TRD-develop
50
6f1e466d 51Revision 1.1.4.1 2000/05/08 15:09:01 cblume
52Introduce AliTRDdigitsManager
53
c0dd96c3 54Revision 1.1 2000/02/28 18:58:54 cblume
55Add new TRD classes
56
f7336fa3 57*/
58
59///////////////////////////////////////////////////////////////////////////////
60// //
61// TRD cluster finder for the slow simulator.
62// //
63///////////////////////////////////////////////////////////////////////////////
64
65#include <TF1.h>
94de3818 66#include <TTree.h>
793ff80c 67#include <TH1.h>
f7336fa3 68
793ff80c 69#include "AliRun.h"
70
71#include "AliTRD.h"
f7336fa3 72#include "AliTRDclusterizerV1.h"
73#include "AliTRDmatrix.h"
74#include "AliTRDgeometry.h"
75#include "AliTRDdigitizer.h"
6f1e466d 76#include "AliTRDdataArrayF.h"
793ff80c 77#include "AliTRDdataArrayI.h"
78#include "AliTRDdigitsManager.h"
f7336fa3 79
80ClassImp(AliTRDclusterizerV1)
81
82//_____________________________________________________________________________
83AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
84{
85 //
86 // AliTRDclusterizerV1 default constructor
87 //
88
6f1e466d 89 fDigitsManager = NULL;
f7336fa3 90
3e1a3ad8 91 fClusMaxThresh = 0;
92 fClusSigThresh = 0;
93
f7336fa3 94}
95
96//_____________________________________________________________________________
97AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
98 :AliTRDclusterizer(name,title)
99{
100 //
101 // AliTRDclusterizerV1 default constructor
102 //
103
6f1e466d 104 fDigitsManager = new AliTRDdigitsManager();
f7336fa3 105
106 Init();
107
108}
109
8230f242 110//_____________________________________________________________________________
dd9a6ee3 111AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
8230f242 112{
113 //
114 // AliTRDclusterizerV1 copy constructor
115 //
116
dd9a6ee3 117 ((AliTRDclusterizerV1 &) c).Copy(*this);
8230f242 118
119}
120
f7336fa3 121//_____________________________________________________________________________
122AliTRDclusterizerV1::~AliTRDclusterizerV1()
123{
8230f242 124 //
125 // AliTRDclusterizerV1 destructor
126 //
f7336fa3 127
6f1e466d 128 if (fDigitsManager) {
129 delete fDigitsManager;
f7336fa3 130 }
131
132}
133
dd9a6ee3 134//_____________________________________________________________________________
135AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
136{
137 //
138 // Assignment operator
139 //
140
141 if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
142 return *this;
143
144}
145
8230f242 146//_____________________________________________________________________________
43da34c0 147void AliTRDclusterizerV1::Copy(TObject &c)
8230f242 148{
149 //
150 // Copy function
151 //
152
43da34c0 153 ((AliTRDclusterizerV1 &) c).fClusMaxThresh = fClusMaxThresh;
154 ((AliTRDclusterizerV1 &) c).fClusSigThresh = fClusSigThresh;
43da34c0 155 ((AliTRDclusterizerV1 &) c).fDigitsManager = NULL;
8230f242 156
157 AliTRDclusterizer::Copy(c);
158
159}
160
f7336fa3 161//_____________________________________________________________________________
162void AliTRDclusterizerV1::Init()
163{
164 //
165 // Initializes the cluster finder
166 //
167
168 // The default parameter for the clustering
3e1a3ad8 169 fClusMaxThresh = 3;
170 fClusSigThresh = 1;
f7336fa3 171
172}
173
174//_____________________________________________________________________________
175Bool_t AliTRDclusterizerV1::ReadDigits()
176{
177 //
178 // Reads the digits arrays from the input aliroot file
179 //
180
181 if (!fInputFile) {
182 printf("AliTRDclusterizerV1::ReadDigits -- ");
183 printf("No input file open\n");
184 return kFALSE;
185 }
186
f7336fa3 187 // Read in the digit arrays
6f1e466d 188 return (fDigitsManager->ReadDigits());
f7336fa3 189
190}
191
192//_____________________________________________________________________________
793ff80c 193Bool_t AliTRDclusterizerV1::MakeClusters()
f7336fa3 194{
195 //
196 // Generates the cluster.
197 //
198
199 Int_t row, col, time;
200
3e1a3ad8 201 if (fTRD->IsVersion() != 1) {
f7336fa3 202 printf("AliTRDclusterizerV1::MakeCluster -- ");
203 printf("TRD must be version 1 (slow simulator).\n");
204 return kFALSE;
205 }
206
207 // Get the geometry
3e1a3ad8 208 AliTRDgeometry *geo = fTRD->GetGeometry();
f7336fa3 209
210 printf("AliTRDclusterizerV1::MakeCluster -- ");
211 printf("Start creating clusters.\n");
212
8230f242 213 AliTRDdataArrayI *digits;
793ff80c 214 AliTRDdataArrayI *track0;
215 AliTRDdataArrayI *track1;
216 AliTRDdataArrayI *track2;
f7336fa3 217
3e1a3ad8 218 // Threshold value for the maximum
219 Int_t maxThresh = fClusMaxThresh;
220 // Threshold value for the digit signal
221 Int_t sigThresh = fClusSigThresh;
f7336fa3 222
223 // Iteration limit for unfolding procedure
8230f242 224 const Float_t kEpsilon = 0.01;
f7336fa3 225
8230f242 226 const Int_t kNclus = 3;
227 const Int_t kNsig = 5;
3e1a3ad8 228 const Int_t kNtrack = 3 * kNclus;
229
230 Float_t padSignal[kNsig];
231 Float_t clusterSignal[kNclus];
232 Float_t clusterPads[kNclus];
233 Int_t clusterDigit[kNclus];
234 Int_t clusterTracks[kNtrack];
f7336fa3 235
236 Int_t chamBeg = 0;
793ff80c 237 Int_t chamEnd = AliTRDgeometry::Ncham();
3e1a3ad8 238 if (fTRD->GetSensChamber() >= 0) {
239 chamBeg = fTRD->GetSensChamber();
6f1e466d 240 chamEnd = chamBeg + 1;
f7336fa3 241 }
242 Int_t planBeg = 0;
793ff80c 243 Int_t planEnd = AliTRDgeometry::Nplan();
3e1a3ad8 244 if (fTRD->GetSensPlane() >= 0) {
245 planBeg = fTRD->GetSensPlane();
f7336fa3 246 planEnd = planBeg + 1;
247 }
248 Int_t sectBeg = 0;
793ff80c 249 Int_t sectEnd = AliTRDgeometry::Nsect();
f7336fa3 250
3e1a3ad8 251 // Start clustering in every chamber
f7336fa3 252 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
253 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
254 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
255
3e1a3ad8 256 if (fTRD->GetSensSector() >= 0) {
257 Int_t sens1 = fTRD->GetSensSector();
258 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
793ff80c 259 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
260 * AliTRDgeometry::Nsect();
dd9a6ee3 261 if (sens1 < sens2) {
9d0b222b 262 if ((isect < sens1) || (isect >= sens2)) continue;
dd9a6ee3 263 }
264 else {
9d0b222b 265 if ((isect < sens1) && (isect >= sens2)) continue;
dd9a6ee3 266 }
9d0b222b 267 }
268
8230f242 269 Int_t idet = geo->GetDetector(iplan,icham,isect);
f7336fa3 270
3e1a3ad8 271 Int_t nClusters = 0;
272 Int_t nClustersUnfold = 0;
273
f7336fa3 274 printf("AliTRDclusterizerV1::MakeCluster -- ");
275 printf("Analyzing chamber %d, plane %d, sector %d.\n"
3e1a3ad8 276 ,icham,iplan,isect);
f7336fa3 277
3e1a3ad8 278 Int_t nRowMax = geo->GetRowMax(iplan,icham,isect);
279 Int_t nColMax = geo->GetColMax(iplan);
280 Int_t nTimeBefore = geo->GetTimeBefore();
281 Int_t nTimeTotal = geo->GetTimeTotal();
f7336fa3 282
3e1a3ad8 283 // Get the digits
8230f242 284 digits = fDigitsManager->GetDigits(idet);
3e1a3ad8 285 digits->Expand();
793ff80c 286 track0 = fDigitsManager->GetDictionary(idet,0);
3e1a3ad8 287 track0->Expand();
793ff80c 288 track1 = fDigitsManager->GetDictionary(idet,1);
3e1a3ad8 289 track1->Expand();
793ff80c 290 track2 = fDigitsManager->GetDictionary(idet,2);
3e1a3ad8 291 track2->Expand();
292
293 // Loop through the chamber and find the maxima
294 for ( row = 0; row < nRowMax; row++) {
295 for ( col = 2; col < nColMax; col++) {
296 for (time = 0; time < nTimeTotal; time++) {
297
298 Int_t signalL = digits->GetDataUnchecked(row,col ,time);
299 Int_t signalM = digits->GetDataUnchecked(row,col-1,time);
300 Int_t signalR = digits->GetDataUnchecked(row,col-2,time);
301
302 // Look for the maximum
303 if ((signalM >= maxThresh) &&
304 (signalL >= sigThresh) &&
305 (signalR >= sigThresh)) {
306 if ((signalL < signalM) && (signalR < signalM)) {
307 // Maximum found, mark the position by a negative signal
308 digits->SetDataUnchecked(row,col-1,time,-signalM);
309 }
310 }
311
312 }
313 }
314 }
315
316 // Now check the maxima and calculate the cluster position
317 for ( row = 0; row < nRowMax ; row++) {
318 for ( col = 1; col < nColMax-1; col++) {
319 for (time = 0; time < nTimeTotal; time++) {
320
321 // Maximum found ?
322 if (digits->GetDataUnchecked(row,col,time) < 0) {
f7336fa3 323
9d0b222b 324 Int_t iPad;
8230f242 325 for (iPad = 0; iPad < kNclus; iPad++) {
3e1a3ad8 326 Int_t iPadCol = col - 1 + iPad;
327 clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
328 ,iPadCol
329 ,time));
330 clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
331 clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
332 clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
333 clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
f7336fa3 334 }
335
3e1a3ad8 336 // Neighbouring maximum on right side?
337 Int_t iType = 0;
338 if (col < nColMax-3) {
339 if (digits->GetDataUnchecked(row,col+2,time) < 0) {
340 for (iPad = 0; iPad < kNsig; iPad++) {
341 padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
342 ,col-1+iPad
343 ,time));
f7336fa3 344 }
3e1a3ad8 345 // Unfold the two maxima and set the signal on
346 // the overlapping pad to the ratio
347 clusterSignal[2] *= Unfold(kEpsilon,padSignal);
348 iType = 1;
f7336fa3 349 }
350 }
f7336fa3 351
3e1a3ad8 352 Float_t clusterCharge = clusterSignal[0]
353 + clusterSignal[1]
354 + clusterSignal[2];
355
356 // Calculate the position of the cluster by using the
357 // center of gravity method
358 clusterPads[0] = row + 0.5;
359 clusterPads[1] = col + 0.5
360 + (clusterSignal[2] - clusterSignal[0])
361 / clusterCharge;
362 // Take the shift of the additional time bins into account
363 clusterPads[2] = time - nTimeBefore + 0.5;
364
365 if (fVerbose) {
366 printf("-----------------------------------------------------------\n");
367 printf("Create cluster no. %d\n",nClusters);
368 printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
369 ,clusterPads[1]
370 ,clusterPads[2]);
371 printf("Indices: %d, %d, %d\n",clusterDigit[0]
372 ,clusterDigit[1]
373 ,clusterDigit[2]);
374 printf("Total charge = %f\n",clusterCharge);
375 printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
376 ,clusterTracks[1]
377 ,clusterTracks[2]);
378 printf(" pad1 %d, %d, %d\n",clusterTracks[3]
379 ,clusterTracks[4]
380 ,clusterTracks[5]);
381 printf(" pad2 %d, %d, %d\n",clusterTracks[6]
382 ,clusterTracks[7]
383 ,clusterTracks[8]);
f7336fa3 384 }
385
3e1a3ad8 386 nClusters++;
387 if (iType == 1) nClustersUnfold++;
f7336fa3 388
389 // Add the cluster to the output array
3e1a3ad8 390 fTRD->AddCluster(clusterPads
391 ,clusterDigit
392 ,idet
393 ,clusterCharge
394 ,clusterTracks
395 ,iType);
f7336fa3 396
397 }
3e1a3ad8 398 }
399 }
400 }
f7336fa3 401
3e1a3ad8 402 // Compress the arrays
403 digits->Compress(1,0);
404 track0->Compress(1,0);
405 track1->Compress(1,0);
406 track2->Compress(1,0);
f7336fa3 407
3e1a3ad8 408 // Write the cluster and reset the array
793ff80c 409 WriteClusters(idet);
3e1a3ad8 410 fTRD->ResetRecPoints();
793ff80c 411
3e1a3ad8 412 printf("AliTRDclusterizerV1::MakeCluster -- ");
413 printf("Found %d (%d unfolded) clusters.\n"
414 ,nClusters,nClustersUnfold);
f7336fa3 415
3e1a3ad8 416 }
417 }
418 }
f7336fa3 419
f7336fa3 420 printf("AliTRDclusterizerV1::MakeCluster -- ");
421 printf("Done.\n");
422
423 return kTRUE;
424
425}
426
427//_____________________________________________________________________________
428Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal)
429{
430 //
431 // Method to unfold neighbouring maxima.
432 // The charge ratio on the overlapping pad is calculated
433 // until there is no more change within the range given by eps.
434 // The resulting ratio is then returned to the calling method.
435 //
436
3e1a3ad8 437 Int_t itStep = 0; // Count iteration steps
f7336fa3 438
3e1a3ad8 439 Float_t ratio = 0.5; // Start value for ratio
440 Float_t prevRatio = 0; // Store previous ratio
f7336fa3 441
3e1a3ad8 442 Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal
443 Float_t newRightSignal[3] = {0}; // Array to store right cluster signal
f7336fa3 444
3e1a3ad8 445 // Start the iteration
f7336fa3 446 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
447
448 itStep++;
449 prevRatio = ratio;
450
3e1a3ad8 451 // Cluster position according to charge ratio
452 Float_t maxLeft = (ratio*padSignal[2] - padSignal[0])
453 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
454 Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
455 / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
f7336fa3 456
3e1a3ad8 457 // Set cluster charge ratio
f7336fa3 458 Float_t ampLeft = padSignal[1];
459 Float_t ampRight = padSignal[3];
460
3e1a3ad8 461 // Apply pad response to parameters
462 newLeftSignal[0] = ampLeft * PadResponse(-1 - maxLeft);
463 newLeftSignal[1] = ampLeft * PadResponse( 0 - maxLeft);
464 newLeftSignal[2] = ampLeft * PadResponse( 1 - maxLeft);
f7336fa3 465
3e1a3ad8 466 newRightSignal[0] = ampRight * PadResponse(-1 - maxRight);
467 newRightSignal[1] = ampRight * PadResponse( 0 - maxRight);
468 newRightSignal[2] = ampRight * PadResponse( 1 - maxRight);
f7336fa3 469
3e1a3ad8 470 // Calculate new overlapping ratio
471 ratio = newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0]);
f7336fa3 472
473 }
474
475 return ratio;
476
477}
478
479//_____________________________________________________________________________
480Float_t AliTRDclusterizerV1::PadResponse(Float_t x)
481{
482 //
483 // The pad response for the chevron pads.
484 // We use a simple Gaussian approximation which should be good
485 // enough for our purpose.
3e1a3ad8 486 // Updated for new PRF 1/5/01.
f7336fa3 487 //
488
489 // The parameters for the response function
3e1a3ad8 490 const Float_t kA = 0.8303;
491 const Float_t kB = -0.00392;
492 const Float_t kC = 0.472 * 0.472;
493 const Float_t kD = 2.19;
f7336fa3 494
3e1a3ad8 495 Float_t pr = kA * (kB + TMath::Exp(-TMath::Power(x*x,kD) / (2.*kC)));
f7336fa3 496
497 return (pr);
498
499}