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