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