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