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