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