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"
34 #include "AliTRDclusterizerMI.h"
35 #include "AliTRDmatrix.h"
36 #include "AliTRDgeometry.h"
37 #include "AliTRDdigitizer.h"
38 #include "AliTRDdataArrayF.h"
39 #include "AliTRDdataArrayI.h"
40 #include "AliTRDdigitsManager.h"
41 #include "AliTRDparameter.h"
42 #include "AliTRDclusterMI.h"
44 ClassImp(AliTRDclusterizerMI)
46 //_____________________________________________________________________________
47 AliTRDclusterizerMI::AliTRDclusterizerMI():AliTRDclusterizerV1()
50 // AliTRDclusterizerMI default constructor
54 //_____________________________________________________________________________
55 AliTRDclusterizerMI::AliTRDclusterizerMI(const Text_t* name, const Text_t* title)
56 :AliTRDclusterizerV1(name,title)
59 // AliTRDclusterizerMI default constructor
63 //_____________________________________________________________________________
64 AliTRDclusterizerMI::~AliTRDclusterizerMI()
67 // AliTRDclusterizerMI destructor
72 AliTRDclusterMI * AliTRDclusterizerMI::AddCluster()
74 AliTRDclusterMI *c = new AliTRDclusterMI();
75 fClusterContainer->Add(c);
79 void AliTRDclusterizerMI::SetCluster(AliTRDclusterMI * cl,Float_t *pos, Int_t det, Float_t amp
80 ,Int_t *tracks, Float_t *sig, Int_t iType, Float_t sigmay, Float_t relpad)
86 cl->AddTrackIndex(tracks);
90 cl->SetSigmaY2(sig[0]);
91 cl->SetSigmaZ2(sig[1]);
92 cl->SetLocalTimeBin(((Int_t) pos[2]));
94 cl->SetRelPos(relpad);
98 void AliTRDclusterizerMI::MakeCluster(Float_t * padSignal, Float_t * pos, Float_t &sigma, Float_t & relpad)
105 Float_t signal[3]={padSignal[0],padSignal[1],padSignal[2]};
107 signal[0] = 0.015*(signal[1]*signal[1])/(signal[2]+0.5);
108 if (signal[0]>2) signal[0]=2;
111 signal[2] = 0.015*(signal[1]*signal[1])/(signal[0]+0.5);
112 if (signal[2]>2) signal[2]=2;
115 for (Int_t i=-1;i<=1;i++){
117 sumx +=signal[i+1]*float(i);
118 sumx2 +=signal[i+1]*float(i)*float(i);
122 sigma = sumx2/sum-pos[0]*pos[0];
126 //_____________________________________________________________________________
127 Bool_t AliTRDclusterizerMI::MakeClusters()
130 // Generates the cluster.
133 //////////////////////
134 //STUPIDITY to be fixed later
135 fClusterContainer = fTRD->RecPoints();
137 Int_t row, col, time;
139 if (fTRD->IsVersion() != 1) {
140 printf("<AliTRDclusterizerMI::MakeCluster> ");
141 printf("TRD must be version 1 (slow simulator).\n");
146 AliTRDgeometry *geo = fTRD->GetGeometry();
148 // Create a default parameter class if none is defined
150 fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter");
151 printf("<AliTRDclusterizerMI::MakeCluster> ");
152 printf("Create the default parameter object.\n");
155 Float_t timeBinSize = fPar->GetTimeBinSize();
156 // Half of ampl.region
157 const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.;
159 Float_t omegaTau = fPar->GetOmegaTau();
161 printf("<AliTRDclusterizerMI::MakeCluster> ");
162 printf("OmegaTau = %f \n",omegaTau);
163 printf("<AliTRDclusterizerMI::MakeCluster> ");
164 printf("Start creating clusters.\n");
167 AliTRDdataArrayI *digits;
168 AliTRDdataArrayI *track0;
169 AliTRDdataArrayI *track1;
170 AliTRDdataArrayI *track2;
172 // Threshold value for the maximum
173 Int_t maxThresh = fPar->GetClusMaxThresh();
174 // Threshold value for the digit signal
175 Int_t sigThresh = fPar->GetClusSigThresh();
177 // Iteration limit for unfolding procedure
178 const Float_t kEpsilon = 0.01;
180 const Int_t kNclus = 3;
181 const Int_t kNsig = 5;
182 const Int_t kNtrack = 3 * kNclus;
187 Float_t ratioLeft = 1.0;
188 Float_t ratioRight = 1.0;
190 Float_t padSignal[kNsig];
191 Float_t clusterSignal[kNclus];
192 Float_t clusterPads[kNclus];
193 Int_t clusterDigit[kNclus];
194 Int_t clusterTracks[kNtrack];
197 Int_t chamEnd = AliTRDgeometry::Ncham();
198 if (fTRD->GetSensChamber() >= 0) {
199 chamBeg = fTRD->GetSensChamber();
200 chamEnd = chamBeg + 1;
203 Int_t planEnd = AliTRDgeometry::Nplan();
204 if (fTRD->GetSensPlane() >= 0) {
205 planBeg = fTRD->GetSensPlane();
206 planEnd = planBeg + 1;
209 Int_t sectEnd = AliTRDgeometry::Nsect();
211 // Start clustering in every chamber
212 for (Int_t icham = chamBeg; icham < chamEnd; icham++) {
213 for (Int_t iplan = planBeg; iplan < planEnd; iplan++) {
214 for (Int_t isect = sectBeg; isect < sectEnd; isect++) {
216 if (fTRD->GetSensSector() >= 0) {
217 Int_t sens1 = fTRD->GetSensSector();
218 Int_t sens2 = sens1 + fTRD->GetSensSectorRange();
219 sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect()))
220 * AliTRDgeometry::Nsect();
222 if ((isect < sens1) || (isect >= sens2)) continue;
225 if ((isect < sens1) && (isect >= sens2)) continue;
229 Int_t idet = geo->GetDetector(iplan,icham,isect);
232 Int_t nClusters2pad = 0;
233 Int_t nClusters3pad = 0;
234 Int_t nClusters4pad = 0;
235 Int_t nClusters5pad = 0;
236 Int_t nClustersLarge = 0;
239 printf("<AliTRDclusterizerMI::MakeCluster> ");
240 printf("Analyzing chamber %d, plane %d, sector %d.\n"
244 Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect);
245 Int_t nColMax = fPar->GetColMax(iplan);
246 Int_t nTimeBefore = fPar->GetTimeBefore();
247 Int_t nTimeTotal = fPar->GetTimeTotal();
249 Float_t row0 = fPar->GetRow0(iplan,icham,isect);
250 Float_t col0 = fPar->GetCol0(iplan);
251 Float_t rowSize = fPar->GetRowPadSize(iplan,icham,isect);
252 Float_t colSize = fPar->GetColPadSize(iplan);
255 digits = fDigitsManager->GetDigits(idet);
257 track0 = fDigitsManager->GetDictionary(idet,0);
259 track1 = fDigitsManager->GetDictionary(idet,1);
261 track2 = fDigitsManager->GetDictionary(idet,2);
264 // Loop through the chamber and find the maxima
265 for ( row = 0; row < nRowMax; row++) {
266 for ( col = 2; col < nColMax; col++) {
267 for (time = 0; time < nTimeTotal; time++) {
269 Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time));
270 Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
271 Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
273 // Look for the maximum
274 if (signalM >= maxThresh) {
275 if (((signalL >= sigThresh) &&
276 (signalL < signalM)) ||
277 ((signalR >= sigThresh) &&
278 (signalR < signalM))) {
279 // Maximum found, mark the position by a negative signal
280 digits->SetDataUnchecked(row,col-1,time,-signalM);
288 // Now check the maxima and calculate the cluster position
289 for ( row = 0; row < nRowMax ; row++) {
290 for (time = 0; time < nTimeTotal; time++) {
291 for ( col = 1; col < nColMax-1; col++) {
294 if (digits->GetDataUnchecked(row,col,time) < 0) {
297 for (iPad = 0; iPad < kNclus; iPad++) {
298 Int_t iPadCol = col - 1 + iPad;
299 clusterSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
302 clusterDigit[iPad] = digits->GetIndexUnchecked(row,iPadCol,time);
303 clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1;
304 clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1;
305 clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1;
308 // Count the number of pads in the cluster
311 while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time))
315 if (col-ii < 0) break;
318 while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time))
322 if (col+ii+1 >= nColMax) break;
349 // Don't analyze large clusters
350 //if (iType == 4) continue;
352 // Look for 5 pad cluster with minimum in the middle
353 Bool_t fivePadCluster = kFALSE;
354 if (col < nColMax-3) {
355 if (digits->GetDataUnchecked(row,col+2,time) < 0) {
356 fivePadCluster = kTRUE;
358 if ((fivePadCluster) && (col < nColMax-5)) {
359 if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) {
360 fivePadCluster = kFALSE;
363 if ((fivePadCluster) && (col > 1)) {
364 if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) {
365 fivePadCluster = kFALSE;
371 // Modify the signal of the overlapping pad for the left part
372 // of the cluster which remains from a previous unfolding
374 clusterSignal[0] *= ratioLeft;
379 // Unfold the 5 pad cluster
380 if (fivePadCluster) {
381 for (iPad = 0; iPad < kNsig; iPad++) {
382 padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row
386 // Unfold the two maxima and set the signal on
387 // the overlapping pad to the ratio
388 ratioRight = Unfold(kEpsilon,iplan,padSignal);
389 ratioLeft = 1.0 - ratioRight;
390 clusterSignal[2] *= ratioRight;
395 Float_t clusterCharge = clusterSignal[0]
399 // The position of the cluster
400 clusterPads[0] = row + 0.5;
401 // Take the shift of the additional time bins into account
402 clusterPads[2] = time - nTimeBefore + 0.5;
406 // Calculate the position of the cluster by using the
407 // lookup table method
408 clusterPads[1] = col + 0.5
409 + fPar->LUTposition(iplan,clusterSignal[0]
416 // Calculate the position of the cluster by using the
417 // center of gravity method
418 clusterPads[1] = col + 0.5
419 + (clusterSignal[2] - clusterSignal[0])
424 Float_t q0 = clusterSignal[0];
425 Float_t q1 = clusterSignal[1];
426 Float_t q2 = clusterSignal[2];
427 Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
428 (clusterCharge*clusterCharge);
430 // Correct for ExB displacement
432 Int_t local_time_bin = (Int_t) clusterPads[2];
433 Float_t driftLength = local_time_bin * timeBinSize + kAmWidth;
434 Float_t colSize = fPar->GetColPadSize(iplan);
435 Float_t deltaY = omegaTau*driftLength/colSize;
436 clusterPads[1] = clusterPads[1] - deltaY;
440 // Calculate the position and the error
441 Float_t clusterPos[3];
442 clusterPos[0] = clusterPads[1] * colSize + col0;
443 clusterPos[1] = clusterPads[0] * rowSize + row0;
444 clusterPos[2] = clusterPads[2];
445 Float_t clusterSig[2];
446 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
447 clusterSig[1] = rowSize * rowSize / 12.;
450 AliTRDclusterMI * cluster = AddCluster();
451 Float_t sigma, relpos;
452 MakeCluster(clusterSignal, clusterPos, sigma,relpos);
454 clusterPos[0] = clusterPads[1] * colSize + col0;
455 clusterPos[1] = clusterPads[0] * rowSize + row0;
456 clusterPos[2] = clusterPads[2];
457 SetCluster(cluster, clusterPos,idet,clusterCharge,clusterTracks,clusterSig,iType,sigma,relpos);
458 // Add the cluster to the output array
459 // fTRD->AddCluster(clusterPos
471 // Compress the arrays
472 digits->Compress(1,0);
473 track0->Compress(1,0);
474 track1->Compress(1,0);
475 track2->Compress(1,0);
477 // Write the cluster and reset the array
479 fTRD->ResetRecPoints();
482 printf("<AliTRDclusterizerMI::MakeCluster> ");
483 printf("Found %d clusters in total.\n"
485 printf(" 2pad: %d\n",nClusters2pad);
486 printf(" 3pad: %d\n",nClusters3pad);
487 printf(" 4pad: %d\n",nClusters4pad);
488 printf(" 5pad: %d\n",nClusters5pad);
489 printf(" Large: %d\n",nClustersLarge);
497 printf("<AliTRDclusterizerMI::MakeCluster> ");
505 //_____________________________________________________________________________
506 Float_t AliTRDclusterizerMI::Unfold(Float_t eps, Int_t plane, Float_t* padSignal)
509 // Method to unfold neighbouring maxima.
510 // The charge ratio on the overlapping pad is calculated
511 // until there is no more change within the range given by eps.
512 // The resulting ratio is then returned to the calling method.
516 Int_t itStep = 0; // Count iteration steps
518 Float_t ratio = 0.5; // Start value for ratio
519 Float_t prevRatio = 0; // Store previous ratio
521 Float_t newLeftSignal[3] = {0}; // Array to store left cluster signal
522 Float_t newRightSignal[3] = {0}; // Array to store right cluster signal
523 Float_t newSignal[3] = {0};
525 // Start the iteration
526 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
531 // Cluster position according to charge ratio
532 Float_t maxLeft = (ratio*padSignal[2] - padSignal[0])
533 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
534 Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
535 / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
537 // Set cluster charge ratio
538 irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal);
539 Float_t ampLeft = padSignal[1] / newSignal[1];
540 irc = fPar->PadResponse(1.0,maxRight,plane,newSignal);
541 Float_t ampRight = padSignal[3] / newSignal[1];
543 // Apply pad response to parameters
544 irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
545 irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal);
547 // Calculate new overlapping ratio
548 ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] /
549 (newLeftSignal[2] + newRightSignal[0]));