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