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 ///////////////////////////////////////////////////////////////////////////////
20 // TRD cluster finder for the slow simulator.
22 ///////////////////////////////////////////////////////////////////////////////
30 #include "AliRunLoader.h"
31 #include "AliLoader.h"
33 #include "AliTRDclusterizerMI.h"
34 #include "AliTRDmatrix.h"
35 #include "AliTRDgeometry.h"
36 #include "AliTRDdataArrayF.h"
37 #include "AliTRDdataArrayI.h"
38 #include "AliTRDdigitsManager.h"
39 #include "AliTRDparameter.h"
40 #include "AliTRDclusterMI.h"
41 #include "AliTRDpadPlane.h"
43 ClassImp(AliTRDclusterizerMI)
45 //_____________________________________________________________________________
46 AliTRDclusterizerMI::AliTRDclusterizerMI():AliTRDclusterizerV1()
49 // AliTRDclusterizerMI default constructor
53 //_____________________________________________________________________________
54 AliTRDclusterizerMI::AliTRDclusterizerMI(const Text_t* name, const Text_t* title)
55 :AliTRDclusterizerV1(name,title)
58 // AliTRDclusterizerMI default constructor
62 //_____________________________________________________________________________
63 AliTRDclusterizerMI::~AliTRDclusterizerMI()
66 // AliTRDclusterizerMI destructor
70 //_____________________________________________________________________________
71 AliTRDclusterMI * AliTRDclusterizerMI::AddCluster()
77 AliTRDclusterMI *c = new AliTRDclusterMI();
78 fClusterContainer->Add(c);
83 //_____________________________________________________________________________
84 void AliTRDclusterizerMI::SetCluster(AliTRDclusterMI * cl,Double_t *pos
85 , Int_t det, Double_t amp
86 , Int_t *tracks, Double_t *sig, Int_t iType
87 , Double_t sigmay, Double_t relpad)
94 cl->AddTrackIndex(tracks);
98 cl->SetSigmaY2(sig[0]);
99 cl->SetSigmaZ2(sig[1]);
100 cl->SetLocalTimeBin(((Int_t) pos[2]));
102 cl->SetRelPos(relpad);
106 //_____________________________________________________________________________
107 void AliTRDclusterizerMI::MakeCluster(Double_t * padSignal, Double_t * pos
108 , Double_t &sigma, Double_t & relpad)
111 // Does something with the cluster
117 Double_t signal[3]={padSignal[0],padSignal[1],padSignal[2]};
119 signal[0] = 0.015*(signal[1]*signal[1])/(signal[2]+0.5);
120 if (signal[0]>2) signal[0]=2;
123 signal[2] = 0.015*(signal[1]*signal[1])/(signal[0]+0.5);
124 if (signal[2]>2) signal[2]=2;
127 for (Int_t i=-1;i<=1;i++){
129 sumx +=signal[i+1]*float(i);
130 sumx2 +=signal[i+1]*float(i)*float(i);
134 sigma = sumx2/sum-pos[0]*pos[0];
138 //_____________________________________________________________________________
139 Bool_t AliTRDclusterizerMI::MakeClusters()
142 // Generates the cluster.
145 //////////////////////
146 //STUPIDITY to be fixed later
147 fClusterContainer = RecPoints();
149 Int_t row, col, time;
152 if (fTRD->IsVersion() != 1) {
153 printf("<AliTRDclusterizerMI::MakeCluster> ");
154 printf("TRD must be version 1 (slow simulator).\n");
160 AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
162 // Create a default parameter class if none is defined
164 fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
165 printf("<AliTRDclusterizerMI::MakeCluster> ");
166 printf("Create the default parameter object.\n");
169 Double_t timeBinSize = fPar->GetDriftVelocity()
170 / fPar->GetSamplingFrequency();
171 // Half of ampl.region
172 const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
174 Float_t omegaTau = fPar->GetOmegaTau();
176 printf("<AliTRDclusterizerMI::MakeCluster> ");
177 printf("OmegaTau = %f \n",omegaTau);
178 printf("<AliTRDclusterizerMI::MakeCluster> ");
179 printf("Start creating clusters.\n");
182 AliTRDdataArrayI *digits;
183 AliTRDdataArrayI *track0;
184 AliTRDdataArrayI *track1;
185 AliTRDdataArrayI *track2;
187 // Threshold value for the maximum
188 Int_t maxThresh = fPar->GetClusMaxThresh();
189 // Threshold value for the digit signal
190 Int_t sigThresh = fPar->GetClusSigThresh();
192 // Iteration limit for unfolding procedure
193 const Double_t kEpsilon = 0.01;
195 const Int_t kNclus = 3;
196 const Int_t kNsig = 5;
197 const Int_t kNtrack = 3 * kNclus;
202 Double_t ratioLeft = 1.0;
203 Double_t ratioRight = 1.0;
205 Double_t padSignal[kNsig];
206 Double_t clusterSignal[kNclus];
207 Double_t clusterPads[kNclus];
208 Int_t clusterDigit[kNclus];
209 Int_t clusterTracks[kNtrack];
212 Int_t chamEnd = AliTRDgeometry::Ncham();
214 Int_t planEnd = AliTRDgeometry::Nplan();
216 Int_t sectEnd = AliTRDgeometry::Nsect();
218 // Start clustering in every chamber
219 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
220 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
221 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
223 Int_t idet = geo->GetDetector(iplan,icham,isect);
226 Int_t nClusters2pad = 0;
227 Int_t nClusters3pad = 0;
228 Int_t nClusters4pad = 0;
229 Int_t nClusters5pad = 0;
230 Int_t nClustersLarge = 0;
233 printf("<AliTRDclusterizerMI::MakeCluster> ");
234 printf("Analyzing chamber %d, plane %d, sector %d.\n"
238 Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect);
239 Int_t nColMax = fPar->GetColMax(iplan);
240 Int_t nTimeBefore = fPar->GetTimeBefore();
241 Int_t nTimeTotal = fPar->GetTimeTotal();
243 AliTRDpadPlane *padPlane = fPar->GetPadPlane(iplan,icham);
246 digits = fDigitsManager->GetDigits(idet);
248 track0 = fDigitsManager->GetDictionary(idet,0);
250 track1 = fDigitsManager->GetDictionary(idet,1);
252 track2 = fDigitsManager->GetDictionary(idet,2);
255 // Loop through the chamber and find the maxima
256 for ( row = 0; row < nRowMax; row++) {
257 for ( col = 2; col < nColMax; col++) {
258 for (time = 0; time < nTimeTotal; time++) {
260 Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time));
261 Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
262 Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
264 // Look for the maximum
265 if (signalM >= maxThresh) {
266 if (((signalL >= sigThresh) &&
267 (signalL < signalM)) ||
268 ((signalR >= sigThresh) &&
269 (signalR < signalM))) {
270 // Maximum found, mark the position by a negative signal
271 digits->SetDataUnchecked(row,col-1,time,-signalM);
279 // Now check the maxima and calculate the cluster position
280 for ( row = 0; row < nRowMax ; row++) {
281 for (time = 0; time < nTimeTotal; time++) {
282 for ( col = 1; col < nColMax-1; col++) {
285 if (digits->GetDataUnchecked(row,col,time) < 0) {
288 for (iPad = 0; iPad < kNclus; iPad++) {
289 Int_t iPadCol = col - 1 + iPad;
290 clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
293 clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
294 clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
295 clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
296 clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
299 // Count the number of pads in the cluster
302 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time))
306 if (col-ii < 0) break;
309 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
313 if (col+ii+1 >= nColMax) break;
340 // Don't analyze large clusters
341 //if (iType == 4) continue;
343 // Look for 5 pad cluster with minimum in the middle
344 Bool_t fivePadCluster = kFALSE;
345 if (col < nColMax-3) {
346 if (digits->GetDataUnchecked(row,col+2,time) < 0) {
347 fivePadCluster = kTRUE;
349 if ((fivePadCluster) && (col < nColMax-5)) {
350 if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
351 fivePadCluster = kFALSE;
354 if ((fivePadCluster) && (col > 1)) {
355 if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
356 fivePadCluster = kFALSE;
362 // Modify the signal of the overlapping pad for the left part
363 // of the cluster which remains from a previous unfolding
365 clusterSignal[0] *= ratioLeft;
370 // Unfold the 5 pad cluster
371 if (fivePadCluster) {
372 for (iPad = 0; iPad < kNsig; iPad++) {
373 padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
377 // Unfold the two maxima and set the signal on
378 // the overlapping pad to the ratio
379 ratioRight = Unfold(kEpsilon,iplan,padSignal);
380 ratioLeft = 1.0 - ratioRight;
381 clusterSignal[2] *= ratioRight;
386 Double_t clusterCharge = clusterSignal[0]
390 // The position of the cluster
391 clusterPads[0] = row + 0.5;
392 // Take the shift of the additional time bins into account
393 clusterPads[2] = time - nTimeBefore + 0.5;
397 // Calculate the position of the cluster by using the
398 // lookup table method
399 // clusterPads[1] = col + 0.5
400 // + fPar->LUTposition(iplan,clusterSignal[0]
402 // ,clusterSignal[2]);
404 + fPar->LUTposition(iplan,clusterSignal[0]
411 // Calculate the position of the cluster by using the
412 // center of gravity method
413 // clusterPads[1] = col + 0.5
414 // + (clusterSignal[2] - clusterSignal[0])
417 + (clusterSignal[2] - clusterSignal[0])
422 Double_t q0 = clusterSignal[0];
423 Double_t q1 = clusterSignal[1];
424 Double_t q2 = clusterSignal[2];
425 Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
426 (clusterCharge*clusterCharge);
428 // Calculate the position and the error
429 Double_t clusterPos[3];
430 // clusterPos[0] = clusterPads[1] * colSize + col0;
431 // clusterPos[1] = clusterPads[0] * rowSize + row0;
432 clusterPos[0] = padPlane->GetColPos(col) - clusterPads[1];
433 clusterPos[1] = padPlane->GetRowPos(row) - clusterPads[0];
434 clusterPos[2] = clusterPads[2];
435 Double_t clusterSig[2];
436 Double_t colSize = padPlane->GetColSize(col);
437 Double_t rowSize = padPlane->GetRowSize(row);
438 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
439 clusterSig[1] = rowSize * rowSize / 12.;
441 // Correct for ExB displacement
443 Int_t local_time_bin = (Int_t) clusterPads[2];
444 Double_t driftLength = local_time_bin * timeBinSize + kAmWidth;
445 Double_t deltaY = omegaTau * driftLength;
446 clusterPos[1] = clusterPos[1] - deltaY;
451 AliTRDclusterMI * cluster = AddCluster();
452 Double_t sigma, relpos;
453 MakeCluster(clusterSignal, clusterPos, sigma,relpos);
455 // clusterPos[0] = clusterPads[1] * colSize + col0;
456 // clusterPos[1] = clusterPads[0] * rowSize + row0;
457 clusterPos[2] = clusterPads[2];
458 clusterPos[0] = padPlane->GetColPos(col) - clusterPads[1];
459 clusterPos[1] = padPlane->GetRowPos(row) - clusterPads[0];
460 SetCluster(cluster, clusterPos,idet,clusterCharge,clusterTracks,clusterSig,iType,sigma,relpos);
461 // Add the cluster to the output array
462 // fTRD->AddCluster(clusterPos
474 // Compress the arrays
475 digits->Compress(1,0);
476 track0->Compress(1,0);
477 track1->Compress(1,0);
478 track2->Compress(1,0);
480 // Write the cluster and reset the array
485 printf("<AliTRDclusterizerMI::MakeCluster> ");
486 printf("Found %d clusters in total.\n"
488 printf(" 2pad: %d\n",nClusters2pad);
489 printf(" 3pad: %d\n",nClusters3pad);
490 printf(" 4pad: %d\n",nClusters4pad);
491 printf(" 5pad: %d\n",nClusters5pad);
492 printf(" Large: %d\n",nClustersLarge);
500 printf("<AliTRDclusterizerMI::MakeCluster> ");
508 //_____________________________________________________________________________
509 Double_t AliTRDclusterizerMI::Unfold(Double_t eps, Int_t plane, Double_t* padSignal)
512 // Method to unfold neighbouring maxima.
513 // The charge ratio on the overlapping pad is calculated
514 // until there is no more change within the range given by eps.
515 // The resulting ratio is then returned to the calling method.
519 Int_t itStep = 0; // Count iteration steps
521 Double_t ratio = 0.5; // Start value for ratio
522 Double_t prevRatio = 0; // Store previous ratio
524 Double_t newLeftSignal[3] = {0}; // Array to store left cluster signal
525 Double_t newRightSignal[3] = {0}; // Array to store right cluster signal
526 Double_t newSignal[3] = {0};
528 // Start the iteration
529 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
534 // Cluster position according to charge ratio
535 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
536 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
537 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
538 / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
540 // Set cluster charge ratio
541 irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
542 Double_t ampLeft = padSignal[1] / newSignal[1];
543 irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
544 Double_t ampRight = padSignal[3] / newSignal[1];
546 // Apply pad response to parameters
547 irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
548 irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
550 // Calculate new overlapping ratio
551 ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] /
552 (newLeftSignal[2] + newRightSignal[0]));