Correct response function in the clusterization (M.Ivanov)
[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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 // TRD cluster finder for the slow simulator. 
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include <TF1.h>
25 #include <TTree.h>
26 #include <TH1.h>
27 #include <TFile.h>
28
29 #include "AliRun.h"
30 #include "AliRunLoader.h"
31 #include "AliLoader.h"
32 #include "AliRawReader.h"
33
34 #include "AliTRDclusterizerV1.h"
35 #include "AliTRDmatrix.h"
36 #include "AliTRDgeometry.h"
37 #include "AliTRDdataArrayF.h"
38 #include "AliTRDdataArrayI.h"
39 #include "AliTRDdigitsManager.h"
40 #include "AliTRDparameter.h"
41 #include "AliTRDpadPlane.h"
42 #include "AliTRDrawData.h"
43 #include "AliTRDcalibDB.h"
44 #include "AliTRDRecParam.h"
45 #include "AliTRDCommonParam.h"
46 #include "AliTRDcluster.h"
47
48 ClassImp(AliTRDclusterizerV1)
49
50 //_____________________________________________________________________________
51 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
52 {
53   //
54   // AliTRDclusterizerV1 default constructor
55   //
56
57   fDigitsManager = 0;
58
59 }
60
61 //_____________________________________________________________________________
62 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
63                     :AliTRDclusterizer(name,title)
64 {
65   //
66   // AliTRDclusterizerV1 default constructor
67   //
68
69   fDigitsManager = new AliTRDdigitsManager();
70   fDigitsManager->CreateArrays();
71
72 }
73
74 //_____________________________________________________________________________
75 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
76 :AliTRDclusterizer(c)
77 {
78   //
79   // AliTRDclusterizerV1 copy constructor
80   //
81
82   ((AliTRDclusterizerV1 &) c).Copy(*this);
83
84 }
85
86 //_____________________________________________________________________________
87 AliTRDclusterizerV1::~AliTRDclusterizerV1()
88 {
89   //
90   // AliTRDclusterizerV1 destructor
91   //
92
93   if (fDigitsManager) {
94     delete fDigitsManager;
95     fDigitsManager = NULL;
96   }
97
98 }
99
100 //_____________________________________________________________________________
101 AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
102 {
103   //
104   // Assignment operator
105   //
106
107   if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
108   return *this;
109
110 }
111
112 //_____________________________________________________________________________
113 void AliTRDclusterizerV1::Copy(TObject &c) const
114 {
115   //
116   // Copy function
117   //
118
119   ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
120
121   AliTRDclusterizer::Copy(c);
122
123 }
124
125 //_____________________________________________________________________________
126 Bool_t AliTRDclusterizerV1::ReadDigits()
127 {
128   //
129   // Reads the digits arrays from the input aliroot file
130   //
131
132   if (!fRunLoader) {
133     printf("<AliTRDclusterizerV1::ReadDigits> ");
134     printf("No input file open\n");
135     return kFALSE;
136   }
137   AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
138   if (!loader->TreeD()) loader->LoadDigits();
139
140   // Read in the digit arrays
141   return (fDigitsManager->ReadDigits(loader->TreeD()));
142
143 }
144
145 //_____________________________________________________________________________
146 Bool_t AliTRDclusterizerV1::ReadDigits(AliRawReader* rawReader)
147 {
148   //
149   // Reads the digits arrays from the ddl file
150   //
151
152   AliTRDrawData *raw = new AliTRDrawData();
153   raw->SetDebug(1);
154
155   fDigitsManager = raw->Raw2Digits(rawReader);
156
157   return kTRUE;
158
159 }
160
161 //_____________________________________________________________________________
162 Bool_t AliTRDclusterizerV1::MakeClusters()
163 {
164   //
165   // Generates the cluster.
166   //
167
168   Int_t row, col, time;
169
170   /*
171   if (fTRD->IsVersion() != 1) {
172     printf("<AliTRDclusterizerV1::MakeCluster> ");
173     printf("TRD must be version 1 (slow simulator).\n");
174     return kFALSE; 
175   }
176   */
177
178   // Get the geometry
179   AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
180
181   // Create a default parameter class if none is defined
182   if (!fPar) {
183     fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
184     printf("<AliTRDclusterizerV1::MakeCluster> ");
185     printf("Create the default parameter object.\n");
186   }
187   fPar->Init();
188   
189   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
190   if (!calibration)
191   {
192     printf("<AliTRDclusterizerMI::MakeCluster> ");
193     printf("ERROR getting instance of AliTRDcalibDB");
194     return kFALSE;  
195   }
196   
197   AliTRDRecParam* recParam = AliTRDRecParam::Instance();
198   if (!recParam)
199   {
200     printf("<AliTRDclusterizerMI::MakeCluster> ");
201     printf("ERROR getting instance of AliTRDRecParam");
202     return kFALSE;  
203   }
204   
205   AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance();
206   if (!commonParam)
207   {
208     printf("<AliTRDdigitizer::MakeDigits> ");
209     printf("Could not get common params\n");
210     return kFALSE;
211   }
212     
213   //Float_t timeBinSize = fPar->GetDriftVelocity()
214   //                    / fPar->GetSamplingFrequency();
215   // Half of ampl.region
216   //  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
217
218   //Float_t omegaTau = fPar->GetOmegaTau();
219   if (fVerbose > 0) {
220     //printf("<AliTRDclusterizerV1::MakeCluster> ");
221     //printf("OmegaTau = %f \n",omegaTau);
222     printf("<AliTRDclusterizerV1::MakeCluster> ");
223     printf("Start creating clusters.\n");
224   } 
225
226   AliTRDdataArrayI *digits;
227   AliTRDdataArrayI *track0;
228   AliTRDdataArrayI *track1;
229   AliTRDdataArrayI *track2; 
230
231   // Threshold value for the maximum
232   Int_t maxThresh = recParam->GetClusMaxThresh();   
233   // Threshold value for the digit signal
234   Int_t sigThresh = recParam->GetClusSigThresh();   
235   // Iteration limit for unfolding procedure
236   const Float_t kEpsilon = 0.01;             
237
238   const Int_t   kNclus   = 3;  
239   const Int_t   kNsig    = 5;
240   const Int_t   kNtrack  = 3 * kNclus;
241
242   Int_t    iType         = 0;
243   Int_t    iUnfold       = 0;  
244   Double_t ratioLeft     = 1.0;
245   Double_t ratioRight    = 1.0;
246
247   //
248   Double_t padSignal[kNsig];   
249   Double_t clusterSignal[kNclus];
250   Double_t clusterPads[kNclus];   
251   Int_t    clusterDigit[kNclus];
252   Int_t    clusterTracks[kNtrack];   
253
254   Int_t    chamBeg = 0;
255   Int_t    chamEnd = AliTRDgeometry::Ncham();
256   Int_t    planBeg = 0;
257   Int_t    planEnd = AliTRDgeometry::Nplan();
258   Int_t    sectBeg = 0;
259   Int_t    sectEnd = AliTRDgeometry::Nsect();
260
261   // Start clustering in every chamber
262   for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
263     for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
264       for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
265
266         Int_t idet = geo->GetDetector(iplan,icham,isect);
267
268         Int_t nClusters      = 0;
269         Int_t nClusters2pad  = 0;
270         Int_t nClusters3pad  = 0;
271         Int_t nClusters4pad  = 0;
272         Int_t nClusters5pad  = 0;
273         Int_t nClustersLarge = 0;
274
275         if (fVerbose > 0) {
276           printf("<AliTRDclusterizerV1::MakeCluster> ");
277           printf("Analyzing chamber %d, plane %d, sector %d.\n"
278                 ,icham,iplan,isect);
279         }
280
281         Int_t    nRowMax     = commonParam->GetRowMax(iplan,icham,isect);
282         Int_t    nColMax     = commonParam->GetColMax(iplan);
283         Int_t    nTimeTotal  = calibration->GetNumberOfTimeBins();
284
285         AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham);
286
287         // Get the digits
288         digits = fDigitsManager->GetDigits(idet);
289         digits->Expand();
290         track0 = fDigitsManager->GetDictionary(idet,0);
291         track0->Expand();
292         track1 = fDigitsManager->GetDictionary(idet,1);
293         track1->Expand();
294         track2 = fDigitsManager->GetDictionary(idet,2); 
295         track2->Expand();
296
297         // Loop through the chamber and find the maxima 
298         for ( row = 0;  row <  nRowMax;    row++) {
299           for ( col = 2;  col <  nColMax;    col++) {
300             //for ( col = 4;  col <  nColMax-2;    col++) {
301             for (time = 0; time < nTimeTotal; time++) {
302
303               Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col  ,time));
304               Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
305               Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
306  
307 //            // Look for the maximum
308 //               if (signalM >= maxThresh) {
309 //                 if (((signalL >= sigThresh) &&
310 //                      (signalL <  signalM))  ||
311 //                     ((signalR >= sigThresh) &&
312 //                      (signalR <  signalM))) {
313 //                   // Maximum found, mark the position by a negative signal
314 //                   digits->SetDataUnchecked(row,col-1,time,-signalM);
315 //              }
316 //            }
317               // Look for the maximum
318               if (signalM >= maxThresh) {
319                 if ( (TMath::Abs(signalL)<=signalM) && (TMath::Abs(signalR)<=signalM) && 
320                      (TMath::Abs(signalL)+TMath::Abs(signalR))>sigThresh ) {
321                   // Maximum found, mark the position by a negative signal
322                   digits->SetDataUnchecked(row,col-1,time,-signalM);
323                 }
324               }
325
326             }  
327           }    
328         }      
329
330         // Now check the maxima and calculate the cluster position
331         for ( row = 0;  row <  nRowMax  ;  row++) {
332           for (time = 0; time < nTimeTotal; time++) {
333             for ( col = 1;  col <  nColMax-1;  col++) {
334
335               // Maximum found ?             
336               if (digits->GetDataUnchecked(row,col,time) < 0) {
337
338                 Int_t iPad;
339                 for (iPad = 0; iPad < kNclus; iPad++) {
340                   Int_t iPadCol = col - 1 + iPad;
341                   clusterSignal[iPad]     = TMath::Abs(digits->GetDataUnchecked(row
342                                                                                ,iPadCol
343                                                                                ,time));
344                   clusterDigit[iPad]      = digits->GetIndexUnchecked(row,iPadCol,time);
345                   clusterTracks[3*iPad  ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
346                   clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
347                   clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
348                 }
349
350                 // Count the number of pads in the cluster
351                 Int_t nPadCount = 0;
352                 Int_t ii        = 0;
353                 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii  ,time))
354                                                                   >= sigThresh) {
355                   nPadCount++;
356                   ii++;
357                   if (col-ii   <        0) break;
358                 }
359                 ii = 0;
360                 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
361                                                                   >= sigThresh) {
362                   nPadCount++;
363                   ii++;
364                   if (col+ii+1 >= nColMax) break;
365                 }
366
367                 nClusters++;
368                 switch (nPadCount) {
369                 case 2:
370                   iType = 0;
371                   nClusters2pad++;
372                   break;
373                 case 3:
374                   iType = 1;
375                   nClusters3pad++;
376                   break;
377                 case 4:
378                   iType = 2;
379                   nClusters4pad++;
380                   break;
381                 case 5:
382                   iType = 3;
383                   nClusters5pad++;
384                   break;
385                 default:
386                   iType = 4;
387                   nClustersLarge++;
388                   break;
389                 };
390
391                  // Look for 5 pad cluster with minimum in the middle
392                 Bool_t fivePadCluster = kFALSE;
393                 if (col < nColMax-3) {
394                   if (digits->GetDataUnchecked(row,col+2,time) < 0) {
395                     fivePadCluster = kTRUE;
396                   }
397                   if ((fivePadCluster) && (col < nColMax-5)) {
398                     if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
399                       fivePadCluster = kFALSE;
400                     }
401                   }
402                   if ((fivePadCluster) && (col >         1)) {
403                     if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
404                       fivePadCluster = kFALSE;
405                     }
406                   }
407                 }
408
409                 // 5 pad cluster
410                 // Modify the signal of the overlapping pad for the left part 
411                 // of the cluster which remains from a previous unfolding
412                 if (iUnfold) {
413                   clusterSignal[0] *= ratioLeft;
414                   iType   = 5;
415                   iUnfold = 0;
416                 }
417
418                 // Unfold the 5 pad cluster
419                 if (fivePadCluster) {
420                   for (iPad = 0; iPad < kNsig; iPad++) {
421                     padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
422                                                                          ,col-1+iPad
423                                                                          ,time));
424                   }
425                   // Unfold the two maxima and set the signal on 
426                   // the overlapping pad to the ratio
427                   ratioRight        = Unfold(kEpsilon,iplan,padSignal);
428                   ratioLeft         = 1.0 - ratioRight; 
429                   clusterSignal[2] *= ratioRight;
430                   iType   = 5;
431                   iUnfold = 1;
432                 }
433
434
435                 Double_t clusterCharge = clusterSignal[0]
436                                        + clusterSignal[1]
437                                        + clusterSignal[2];
438                 
439                 // The position of the cluster
440                 clusterPads[0] = row + 0.5;
441                 // Take the shift of the additional time bins into account
442                 clusterPads[2] = time + 0.5;
443
444                 // correct for t0
445                 clusterPads[2] -= calibration->GetT0(idet, col, row);
446                 
447                 if (recParam->LUTOn()) {
448                   // Calculate the position of the cluster by using the
449                   // lookup table method
450                   clusterPads[1] =
451                       recParam->LUTposition(iplan,clusterSignal[0]
452                                                           ,clusterSignal[1]
453                                                           ,clusterSignal[2]);
454                 }
455                 else {
456                   // Calculate the position of the cluster by using the
457                   // center of gravity method
458                   for (Int_t i=0;i<5;i++) padSignal[i]=0;
459                   padSignal[2] = TMath::Abs(digits->GetDataUnchecked(row,col,time));   // central  pad
460                   padSignal[1] = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); // left     pad
461                   padSignal[3] = TMath::Abs(digits->GetDataUnchecked(row,col+1,time)); // right    pad
462                   if (col>2 &&TMath::Abs(digits->GetDataUnchecked(row,col-2,time)<padSignal[1])){
463                     padSignal[0] = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
464                   }
465                   if (col<nColMax-3 &&TMath::Abs(digits->GetDataUnchecked(row,col+2,time)<padSignal[3])){
466                     padSignal[4] = TMath::Abs(digits->GetDataUnchecked(row,col+2,time));
467                   }               
468                   clusterPads[1] =  GetCOG(padSignal);
469                   //Double_t check = fPar->LUTposition(iplan,clusterSignal[0]
470                   //                                        ,clusterSignal[1]
471                   //                                        ,clusterSignal[2]);
472                   //              Float_t diff = clusterPads[1] -  check;
473
474                 }
475
476                 Double_t q0 = clusterSignal[0];
477                 Double_t q1 = clusterSignal[1];
478                 Double_t q2 = clusterSignal[2];
479                 Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
480                                           (clusterCharge*clusterCharge);
481
482                 Float_t vdrift = calibration->GetVdrift(idet, col, row);  
483                 
484                 // Calculate the position and the error
485                 Double_t colSize = padPlane->GetColSize(col);
486                 Double_t rowSize = padPlane->GetRowSize(row);
487                 Double_t clusterPos[3];
488                 clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1]+0.5)*colSize;  // MI change
489                 clusterPos[1] = padPlane->GetRowPos(row) -0.5*rowSize; //MI change
490                 clusterPos[2] = CalcXposFromTimebin(clusterPads[2], vdrift);
491                 Double_t clusterSig[2];
492                 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
493                 clusterSig[1] = rowSize * rowSize / 12.;                                       
494                 
495                 
496                 // Add the cluster to the output array 
497                 AliTRDcluster * cluster = AddCluster(clusterPos
498                           ,(Int_t) clusterPads[2]
499                           ,idet
500                           ,clusterCharge
501                           ,clusterTracks
502                           ,clusterSig
503                            ,iType,clusterPads[1]);
504                 //
505                 //
506                 Short_t signals[7]={0,0,0,0,0,0,0};
507                 for (Int_t jPad = col-3;jPad<=col+3;jPad++){
508                   if (jPad<0 ||jPad>=nColMax-1) continue;
509                   signals[jPad-col+3] =  TMath::Abs(digits->GetDataUnchecked(row,jPad,time));
510                 }
511                 cluster->SetSignals(signals);
512               }
513             } 
514           }   
515         }     
516
517         // Compress the arrays
518         digits->Compress(1,0);
519         track0->Compress(1,0);
520         track1->Compress(1,0);
521         track2->Compress(1,0);
522
523         // Write the cluster and reset the array
524         WriteClusters(idet);
525         ResetRecPoints();
526       }    
527     }      
528   }        
529
530   if (fVerbose > 0) {
531     printf("<AliTRDclusterizerV1::MakeCluster> ");
532     printf("Done.\n");
533   }
534
535   return kTRUE;
536
537 }
538
539 Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5])
540 {
541   //
542   // get COG position
543   // used for clusters with more than 3 pads - where LUT not applicable
544   Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4];
545   Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum;
546   return res;             
547 }
548
549
550
551 //_____________________________________________________________________________
552 Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal)
553 {
554   //
555   // Method to unfold neighbouring maxima.
556   // The charge ratio on the overlapping pad is calculated
557   // until there is no more change within the range given by eps.
558   // The resulting ratio is then returned to the calling method.
559   //
560
561   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
562   if (!calibration)
563   {
564     printf("<AliTRDclusterizerMI::Unfold> ");
565     printf("ERROR getting instance of AliTRDcalibDB");
566     return kFALSE;  
567   }
568   
569   Int_t   irc                = 0;
570   Int_t   itStep             = 0;      // Count iteration steps
571
572   Double_t ratio             = 0.5;    // Start value for ratio
573   Double_t prevRatio         = 0;      // Store previous ratio
574
575   Double_t newLeftSignal[3]  = {0};    // Array to store left cluster signal
576   Double_t newRightSignal[3] = {0};    // Array to store right cluster signal
577   Double_t newSignal[3]      = {0};
578
579   // Start the iteration
580   while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
581
582     itStep++;
583     prevRatio = ratio;
584
585     // Cluster position according to charge ratio
586     Double_t maxLeft  = (ratio*padSignal[2] - padSignal[0]) 
587                       / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
588     Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) 
589                       / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
590
591     // Set cluster charge ratio
592     irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
593     Double_t ampLeft  = padSignal[1] / newSignal[1];
594     irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
595     Double_t ampRight = padSignal[3] / newSignal[1];
596
597     // Apply pad response to parameters
598     irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
599     irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
600
601     // Calculate new overlapping ratio
602     ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] / 
603                           (newLeftSignal[2] + newRightSignal[0]));
604
605   }
606
607   return ratio;
608
609 }
610