Make code compliant to coding conventions
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
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 Revision 1.4  2000/06/07 16:27:01  cblume
19 Try to remove compiler warnings on Sun and HP
20
21 Revision 1.3  2000/05/08 16:17:27  cblume
22 Merge TRD-develop
23
24 Revision 1.1.4.1  2000/05/08 15:09:01  cblume
25 Introduce AliTRDdigitsManager
26
27 Revision 1.1  2000/02/28 18:58:54  cblume
28 Add new TRD classes
29
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"
45 #include "AliTRDdataArrayF.h"
46
47 ClassImp(AliTRDclusterizerV1)
48
49 //_____________________________________________________________________________
50 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
51 {
52   //
53   // AliTRDclusterizerV1 default constructor
54   //
55
56   fDigitsManager = NULL;
57
58 }
59
60 //_____________________________________________________________________________
61 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
62                     :AliTRDclusterizer(name,title)
63 {
64   //
65   // AliTRDclusterizerV1 default constructor
66   //
67
68   fDigitsManager = new AliTRDdigitsManager();
69
70   Init();
71
72 }
73
74 //_____________________________________________________________________________
75 AliTRDclusterizerV1::AliTRDclusterizerV1(AliTRDclusterizerV1 &c)
76 {
77   //
78   // AliTRDclusterizerV1 copy constructor
79   //
80
81   c.Copy(*this);
82
83 }
84
85 //_____________________________________________________________________________
86 AliTRDclusterizerV1::~AliTRDclusterizerV1()
87 {
88   //
89   // AliTRDclusterizerV1 destructor
90   //
91
92   if (fDigitsManager) {
93     delete fDigitsManager;
94   }
95
96 }
97
98 //_____________________________________________________________________________
99 void 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
114 //_____________________________________________________________________________
115 void 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 //_____________________________________________________________________________
129 Bool_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
141   // Read in the digit arrays
142   return (fDigitsManager->ReadDigits());  
143
144 }
145
146 //_____________________________________________________________________________
147 Bool_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
156   AliTRD *trd = (AliTRD*) gAlice->GetDetector("TRD");
157   if (trd->IsVersion() != 1) {
158     printf("AliTRDclusterizerV1::MakeCluster -- ");
159     printf("TRD must be version 1 (slow simulator).\n");
160     return kFALSE; 
161   }
162
163   // Get the geometry
164   AliTRDgeometry *geo = trd->GetGeometry();
165
166   printf("AliTRDclusterizerV1::MakeCluster -- ");
167   printf("Start creating clusters.\n");
168
169   AliTRDdataArrayI *digits;
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
177   const Float_t kEpsilon = 0.01;             
178
179   const Int_t   kNclus   = 3;  
180   const Int_t   kNsig    = 5;
181
182   Int_t chamBeg = 0;
183   Int_t chamEnd = kNcham;
184   if (trd->GetSensChamber()  >= 0) {
185     chamBeg = trd->GetSensChamber();
186     chamEnd = chamBeg + 1;
187   }
188   Int_t planBeg = 0;
189   Int_t planEnd = kNplan;
190   if (trd->GetSensPlane()    >= 0) {
191     planBeg = trd->GetSensPlane();
192     planEnd = planBeg + 1;
193   }
194   Int_t sectBeg = 0;
195   Int_t sectEnd = kNsect;
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
202         if (trd->GetSensSector() >= 0) {
203           Int_t sens1 = trd->GetSensSector();
204           Int_t sens2 = sens1 + trd->GetSensSectorRange();
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
212         Int_t idet = geo->GetDetector(iplan,icham,isect);
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
219         Int_t   nRowMax  = geo->GetRowMax(iplan,icham,isect);
220         Int_t   nColMax  = geo->GetColMax(iplan);
221         Int_t   nTimeMax = geo->GetTimeMax();
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
231         digits = fDigitsManager->GetDigits(idet);
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
238               Int_t signal = digits->GetData(row,col,time);
239               Int_t index  = digits->GetIndex(row,col,time);
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
284                 Float_t ratio                 =  0;    
285                 // Signals on max and neighbouring pads
286                 Float_t padSignal[kNsig]      = {0};   
287                 // Signals from cluster
288                 Float_t clusterSignal[kNclus] = {0};
289                 // Cluster pad info
290                 Float_t clusterPads[kNclus]   = {0};   
291                 // Cluster digit info
292                 Int_t   clusterDigit[kNclus]  = {0};
293
294                 Int_t iPad;
295                 for (iPad = 0; iPad < kNclus; iPad++) {
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
304                     for (iPad = 0; iPad < 5; iPad++) {
305                       padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time);
306                     }
307
308                     // unfold:
309                     ratio = Unfold(kEpsilon, padSignal);
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]) /
323                                    (clusterSignal[0] + clusterSignal[1] + clusterSignal[2]);
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 
355                 trd->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge);
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"
374         ,trd->RecPoints()->GetEntries());
375
376   // Get the pointer to the cluster branch
377   TTree *clusterTree = gAlice->TreeR(); 
378
379   // Fill the cluster-branch
380   printf("AliTRDclusterizerV1::MakeCluster -- ");
381   printf("Fill the cluster tree.\n");
382   clusterTree->Fill();
383   printf("AliTRDclusterizerV1::MakeCluster -- ");
384   printf("Done.\n");
385
386   return kTRUE;
387
388 }
389
390 //_____________________________________________________________________________
391 Float_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 //_____________________________________________________________________________
443 Float_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
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;
456
457   Float_t pr = kA * (kB + TMath::Exp(-x*x / (2. * kC2)));
458
459   return (pr);
460
461 }