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