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