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