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