1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 Revision 1.15.6.2 2002/07/24 10:09:30 alibrary
21 Revision 1.17 2002/06/12 09:54:35 cblume
22 Update of tracking code provided by Sergei
24 Revision 1.16 2002/03/25 20:01:30 cblume
25 Introduce parameter class
27 Revision 1.15 2001/11/14 12:09:11 cblume
28 Use correct name for digitizer
30 Revision 1.14 2001/11/14 10:50:45 cblume
31 Changes in digits IO. Add merging of summable digits
33 Revision 1.13 2001/05/28 17:07:58 hristov
34 Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh)
36 Revision 1.12 2001/05/21 17:42:58 hristov
37 Constant casted to avoid the ambiguity
39 Revision 1.11 2001/05/21 16:45:47 hristov
40 Last minute changes (C.Blume)
42 Revision 1.10 2001/05/07 08:06:44 cblume
43 Speedup of the code. Create only AliTRDcluster
45 Revision 1.9 2000/11/01 14:53:20 cblume
46 Merge with TRD-develop
48 Revision 1.1.4.5 2000/10/15 23:40:01 cblume
51 Revision 1.1.4.4 2000/10/06 16:49:46 cblume
54 Revision 1.1.4.3 2000/10/04 16:34:58 cblume
55 Replace include files by forward declarations
57 Revision 1.1.4.2 2000/09/22 14:49:49 cblume
58 Adapted to tracking code
60 Revision 1.8 2000/10/02 21:28:19 fca
61 Removal of useless dependecies via forward declarations
63 Revision 1.7 2000/06/27 13:08:50 cblume
64 Changed to Copy(TObject &A) to appease the HP-compiler
66 Revision 1.6 2000/06/09 11:10:07 cblume
67 Compiler warnings and coding conventions, next round
69 Revision 1.5 2000/06/08 18:32:58 cblume
70 Make code compliant to coding conventions
72 Revision 1.4 2000/06/07 16:27:01 cblume
73 Try to remove compiler warnings on Sun and HP
75 Revision 1.3 2000/05/08 16:17:27 cblume
78 Revision 1.1.4.1 2000/05/08 15:09:01 cblume
79 Introduce AliTRDdigitsManager
81 Revision 1.1 2000/02/28 18:58:54 cblume
86 ///////////////////////////////////////////////////////////////////////////////
88 // TRD cluster finder for the slow simulator.
90 ///////////////////////////////////////////////////////////////////////////////
100 #include "AliTRDclusterizerV1.h"
101 #include "AliTRDmatrix.h"
102 #include "AliTRDgeometry.h"
103 #include "AliTRDdigitizer.h"
104 #include "AliTRDdataArrayF.h"
105 #include "AliTRDdataArrayI.h"
106 #include "AliTRDdigitsManager.h"
107 #include "AliTRDparameter.h"
109 ClassImp(AliTRDclusterizerV1)
111 //_____________________________________________________________________________
112 AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer()
115 // AliTRDclusterizerV1 default constructor
122 //_____________________________________________________________________________
123 AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title)
124 :AliTRDclusterizer(name,title)
127 // AliTRDclusterizerV1 default constructor
130 fDigitsManager = new AliTRDdigitsManager();
131 fDigitsManager->CreateArrays();
135 //_____________________________________________________________________________
136 AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
139 // AliTRDclusterizerV1 copy constructor
142 ((AliTRDclusterizerV1 &) c).Copy(*this);
146 //_____________________________________________________________________________
147 AliTRDclusterizerV1::~AliTRDclusterizerV1()
150 // AliTRDclusterizerV1 destructor
153 if (fDigitsManager) {
154 delete fDigitsManager;
155 fDigitsManager = NULL;
160 //_____________________________________________________________________________
161 AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
164 // Assignment operator
167 if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
172 //_____________________________________________________________________________
173 void AliTRDclusterizerV1::Copy(TObject &c)
179 ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
181 AliTRDclusterizer::Copy(c);
185 //_____________________________________________________________________________
186 Bool_t AliTRDclusterizerV1::ReadDigits()
189 // Reads the digits arrays from the input aliroot file
193 printf("<AliTRDclusterizerV1::ReadDigits> ");
194 printf("No input file open\n");
198 fDigitsManager->Open(fInputFile->GetName());
200 // Read in the digit arrays
201 return (fDigitsManager->ReadDigits());
205 //_____________________________________________________________________________
206 Bool_t AliTRDclusterizerV1::MakeClusters()
209 // Generates the cluster.
212 Int_t row, col, time;
214 if (fTRD->IsVersion() != 1) {
215 printf("<AliTRDclusterizerV1::MakeCluster> ");
216 printf("TRD must be version 1 (slow simulator).\n");
221 AliTRDgeometry *geo = fTRD->GetGeometry();
223 // Create a default parameter class if none is defined
225 fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
226 printf("<AliTRDclusterizerV1::MakeCluster> ");
227 printf("Create the default parameter object.\n");
230 Float_t timeBinSize = fPar->GetTimeBinSize();
231 // Half of ampl.region
232 const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
234 Float_t omegaTau = fPar->GetOmegaTau();
236 printf("<AliTRDclusterizerV1::MakeCluster> ");
237 printf("OmegaTau = %f \n",omegaTau);
238 printf("<AliTRDclusterizerV1::MakeCluster> ");
239 printf("Start creating clusters.\n");
242 AliTRDdataArrayI *digits;
243 AliTRDdataArrayI *track0;
244 AliTRDdataArrayI *track1;
245 AliTRDdataArrayI *track2;
247 // Threshold value for the maximum
248 Int_t maxThresh = fPar->GetClusMaxThresh();
249 // Threshold value for the digit signal
250 Int_t sigThresh = fPar->GetClusSigThresh();
252 // Iteration limit for unfolding procedure
253 const Float_t kEpsilon = 0.01;
255 const Int_t kNclus = 3;
256 const Int_t kNsig = 5;
257 const Int_t kNtrack = 3 * kNclus;
262 Float_t ratioLeft = 1.0;
263 Float_t ratioRight = 1.0;
265 Float_t padSignal[kNsig];
266 Float_t clusterSignal[kNclus];
267 Float_t clusterPads[kNclus];
268 Int_t clusterDigit[kNclus];
269 Int_t clusterTracks[kNtrack];
272 Int_t chamEnd = AliTRDgeometry::Ncham();
273 if (fTRD->GetSensChamber() >= 0) {
274 chamBeg = fTRD->GetSensChamber();
275 chamEnd = chamBeg + 1;
278 Int_t planEnd = AliTRDgeometry::Nplan();
279 if (fTRD->GetSensPlane() >= 0) {
280 planBeg = fTRD->GetSensPlane();
281 planEnd = planBeg + 1;
284 Int_t sectEnd = AliTRDgeometry::Nsect();
286 // Start clustering in every chamber
287 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
288 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
289 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
291 if (fTRD->GetSensSector() >= 0) {
292 Int_t sens1 = fTRD->GetSensSector();
293 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
294 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
295 * AliTRDgeometry::Nsect();
297 if ((isect < sens1) || (isect >= sens2)) continue;
300 if ((isect < sens1) && (isect >= sens2)) continue;
304 Int_t idet = geo->GetDetector(iplan,icham,isect);
307 Int_t nClusters2pad = 0;
308 Int_t nClusters3pad = 0;
309 Int_t nClusters4pad = 0;
310 Int_t nClusters5pad = 0;
311 Int_t nClustersLarge = 0;
314 printf("<AliTRDclusterizerV1::MakeCluster> ");
315 printf("Analyzing chamber %d, plane %d, sector %d.\n"
319 Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect);
320 Int_t nColMax = fPar->GetColMax(iplan);
321 Int_t nTimeBefore = fPar->GetTimeBefore();
322 Int_t nTimeTotal = fPar->GetTimeTotal();
324 Float_t row0 = fPar->GetRow0(iplan,icham,isect);
325 Float_t col0 = fPar->GetCol0(iplan);
326 Float_t rowSize = fPar->GetRowPadSize(iplan,icham,isect);
327 Float_t colSize = fPar->GetColPadSize(iplan);
330 digits = fDigitsManager->GetDigits(idet);
332 track0 = fDigitsManager->GetDictionary(idet,0);
334 track1 = fDigitsManager->GetDictionary(idet,1);
336 track2 = fDigitsManager->GetDictionary(idet,2);
339 // Loop through the chamber and find the maxima
340 for ( row = 0; row < nRowMax; row++) {
341 for ( col = 2; col < nColMax; col++) {
342 for (time = 0; time < nTimeTotal; time++) {
344 Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time));
345 Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
346 Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
348 // Look for the maximum
349 if (signalM >= maxThresh) {
350 if (((signalL >= sigThresh) &&
351 (signalL < signalM)) ||
352 ((signalR >= sigThresh) &&
353 (signalR < signalM))) {
354 // Maximum found, mark the position by a negative signal
355 digits->SetDataUnchecked(row,col-1,time,-signalM);
363 // Now check the maxima and calculate the cluster position
364 for ( row = 0; row < nRowMax ; row++) {
365 for (time = 0; time < nTimeTotal; time++) {
366 for ( col = 1; col < nColMax-1; col++) {
369 if (digits->GetDataUnchecked(row,col,time) < 0) {
372 for (iPad = 0; iPad < kNclus; iPad++) {
373 Int_t iPadCol = col - 1 + iPad;
374 clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
377 clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
378 clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
379 clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
380 clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
383 // Count the number of pads in the cluster
386 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time))
390 if (col-ii < 0) break;
393 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
397 if (col+ii+1 >= nColMax) break;
424 // Don't analyze large clusters
425 //if (iType == 4) continue;
427 // Look for 5 pad cluster with minimum in the middle
428 Bool_t fivePadCluster = kFALSE;
429 if (col < nColMax-3) {
430 if (digits->GetDataUnchecked(row,col+2,time) < 0) {
431 fivePadCluster = kTRUE;
433 if ((fivePadCluster) && (col < nColMax-5)) {
434 if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
435 fivePadCluster = kFALSE;
438 if ((fivePadCluster) && (col > 1)) {
439 if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
440 fivePadCluster = kFALSE;
446 // Modify the signal of the overlapping pad for the left part
447 // of the cluster which remains from a previous unfolding
449 clusterSignal[0] *= ratioLeft;
454 // Unfold the 5 pad cluster
455 if (fivePadCluster) {
456 for (iPad = 0; iPad < kNsig; iPad++) {
457 padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
461 // Unfold the two maxima and set the signal on
462 // the overlapping pad to the ratio
463 ratioRight = Unfold(kEpsilon,iplan,padSignal);
464 ratioLeft = 1.0 - ratioRight;
465 clusterSignal[2] *= ratioRight;
470 Float_t clusterCharge = clusterSignal[0]
474 // The position of the cluster
475 clusterPads[0] = row + 0.5;
476 // Take the shift of the additional time bins into account
477 clusterPads[2] = time - nTimeBefore + 0.5;
481 // Calculate the position of the cluster by using the
482 // lookup table method
483 clusterPads[1] = col + 0.5
484 + fPar->LUTposition(iplan,clusterSignal[0]
491 // Calculate the position of the cluster by using the
492 // center of gravity method
493 clusterPads[1] = col + 0.5
494 + (clusterSignal[2] - clusterSignal[0])
499 Float_t q0 = clusterSignal[0];
500 Float_t q1 = clusterSignal[1];
501 Float_t q2 = clusterSignal[2];
502 Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
503 (clusterCharge*clusterCharge);
505 // Correct for ExB displacement
507 Int_t local_time_bin = (Int_t) clusterPads[2];
508 Float_t driftLength = local_time_bin * timeBinSize + kAmWidth;
509 Float_t colSize = fPar->GetColPadSize(iplan);
510 Float_t deltaY = omegaTau*driftLength/colSize;
511 clusterPads[1] = clusterPads[1] - deltaY;
515 printf("-----------------------------------------------------------\n");
516 printf("Create cluster no. %d\n",nClusters);
517 printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
520 printf("Indices: %d, %d, %d\n",clusterDigit[0]
523 printf("Total charge = %f\n",clusterCharge);
524 printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
527 printf(" pad1 %d, %d, %d\n",clusterTracks[3]
530 printf(" pad2 %d, %d, %d\n",clusterTracks[6]
533 printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
536 // Calculate the position and the error
537 Float_t clusterPos[3];
538 clusterPos[0] = clusterPads[1] * colSize + col0;
539 clusterPos[1] = clusterPads[0] * rowSize + row0;
540 clusterPos[2] = clusterPads[2];
541 Float_t clusterSig[2];
542 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
543 clusterSig[1] = rowSize * rowSize / 12.;
545 // Add the cluster to the output array
546 fTRD->AddCluster(clusterPos
558 // Compress the arrays
559 digits->Compress(1,0);
560 track0->Compress(1,0);
561 track1->Compress(1,0);
562 track2->Compress(1,0);
564 // Write the cluster and reset the array
566 fTRD->ResetRecPoints();
569 printf("<AliTRDclusterizerV1::MakeCluster> ");
570 printf("Found %d clusters in total.\n"
572 printf(" 2pad: %d\n",nClusters2pad);
573 printf(" 3pad: %d\n",nClusters3pad);
574 printf(" 4pad: %d\n",nClusters4pad);
575 printf(" 5pad: %d\n",nClusters5pad);
576 printf(" Large: %d\n",nClustersLarge);
584 printf("<AliTRDclusterizerV1::MakeCluster> ");
592 //_____________________________________________________________________________
593 Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, Float_t* padSignal)
596 // Method to unfold neighbouring maxima.
597 // The charge ratio on the overlapping pad is calculated
598 // until there is no more change within the range given by eps.
599 // The resulting ratio is then returned to the calling method.
603 Int_t itStep = 0; // Count iteration steps
605 Float_t ratio = 0.5; // Start value for ratio
606 Float_t prevRatio = 0; // Store previous ratio
608 Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal
609 Float_t newRightSignal[3] = {0}; // Array to store right cluster signal
610 Float_t newSignal[3] = {0};
612 // Start the iteration
613 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
618 // Cluster position according to charge ratio
619 Float_t maxLeft = (ratio*padSignal[2] - padSignal[0])
620 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
621 Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
622 / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
624 // Set cluster charge ratio
625 irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
626 Float_t ampLeft = padSignal[1] / newSignal[1];
627 irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
628 Float_t ampRight = padSignal[3] / newSignal[1];
630 // Apply pad response to parameters
631 irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
632 irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
634 // Calculate new overlapping ratio
635 ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] /
636 (newLeftSignal[2] + newRightSignal[0]));