]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TRD/AliTRDclusterizerV1.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Additional data member to AliCluster requested by Raphaelle
[u/mrichter/AliRoot.git] / TRD / AliTRDclusterizerV1.cxx
CommitLineData
cb86ff6e 1
f7336fa3 2/**************************************************************************
3 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * *
5 * Author: The ALICE Off-line Project. *
6 * Contributors are mentioned in the code where appropriate. *
7 * *
8 * Permission to use, copy, modify and distribute this software and its *
9 * documentation strictly for non-commercial purposes is hereby granted *
10 * without fee, provided that the above copyright notice appears in all *
11 * copies and that both the copyright notice and this permission notice *
12 * appear in the supporting documentation. The authors make no claims *
13 * about the suitability of this software for any purpose. It is *
14 * provided "as is" without express or implied warranty. *
15 **************************************************************************/
16
88cb7938 17/* $Id$ */
f7336fa3 18
19///////////////////////////////////////////////////////////////////////////////
20// //
6d50f529 21// TRD cluster finder //
f7336fa3 22// //
23///////////////////////////////////////////////////////////////////////////////
24
25#include <TF1.h>
94de3818 26#include <TTree.h>
793ff80c 27#include <TH1.h>
a819a5f7 28#include <TFile.h>
f7336fa3 29
88cb7938 30#include "AliRunLoader.h"
31#include "AliLoader.h"
928e9fae 32#include "AliRawReader.h"
6d50f529 33#include "AliLog.h"
793ff80c 34
f7336fa3 35#include "AliTRDclusterizerV1.h"
f7336fa3 36#include "AliTRDgeometry.h"
6f1e466d 37#include "AliTRDdataArrayF.h"
793ff80c 38#include "AliTRDdataArrayI.h"
39#include "AliTRDdigitsManager.h"
a5cadd36 40#include "AliTRDpadPlane.h"
928e9fae 41#include "AliTRDrawData.h"
3551db50 42#include "AliTRDcalibDB.h"
3becff3c 43#include "AliTRDSimParam.h"
3551db50 44#include "AliTRDRecParam.h"
45#include "AliTRDCommonParam.h"
c85a4951 46#include "AliTRDcluster.h"
f7336fa3 47
56178ff4 48#include "Cal/AliTRDCalROC.h"
49#include "Cal/AliTRDCalDet.h"
50
f7336fa3 51ClassImp(AliTRDclusterizerV1)
52
53//_____________________________________________________________________________
6d50f529 54AliTRDclusterizerV1::AliTRDclusterizerV1()
55 :AliTRDclusterizer()
56 ,fDigitsManager(NULL)
f7336fa3 57{
58 //
59 // AliTRDclusterizerV1 default constructor
60 //
61
f7336fa3 62}
63
64//_____________________________________________________________________________
a6dd11e9 65AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t *name, const Text_t *title)
6d50f529 66 :AliTRDclusterizer(name,title)
67 ,fDigitsManager(new AliTRDdigitsManager())
f7336fa3 68{
69 //
6d50f529 70 // AliTRDclusterizerV1 constructor
f7336fa3 71 //
72
17b26de4 73 fDigitsManager->CreateArrays();
f7336fa3 74
75}
76
8230f242 77//_____________________________________________________________________________
dd9a6ee3 78AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c)
6d50f529 79 :AliTRDclusterizer(c)
80 ,fDigitsManager(NULL)
8230f242 81{
82 //
83 // AliTRDclusterizerV1 copy constructor
84 //
85
8230f242 86}
87
f7336fa3 88//_____________________________________________________________________________
89AliTRDclusterizerV1::~AliTRDclusterizerV1()
90{
8230f242 91 //
92 // AliTRDclusterizerV1 destructor
93 //
f7336fa3 94
6f1e466d 95 if (fDigitsManager) {
96 delete fDigitsManager;
abaf1f1d 97 fDigitsManager = NULL;
f7336fa3 98 }
99
100}
101
dd9a6ee3 102//_____________________________________________________________________________
103AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c)
104{
105 //
106 // Assignment operator
107 //
108
109 if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this);
110 return *this;
111
112}
113
8230f242 114//_____________________________________________________________________________
e0d47c25 115void AliTRDclusterizerV1::Copy(TObject &c) const
8230f242 116{
117 //
118 // Copy function
119 //
120
17b26de4 121 ((AliTRDclusterizerV1 &) c).fDigitsManager = 0;
8230f242 122
123 AliTRDclusterizer::Copy(c);
124
125}
126
f7336fa3 127//_____________________________________________________________________________
128Bool_t AliTRDclusterizerV1::ReadDigits()
129{
130 //
131 // Reads the digits arrays from the input aliroot file
132 //
133
88cb7938 134 if (!fRunLoader) {
6d50f529 135 AliError("No run loader available");
f7336fa3 136 return kFALSE;
137 }
6d50f529 138
88cb7938 139 AliLoader* loader = fRunLoader->GetLoader("TRDLoader");
6d50f529 140 if (!loader->TreeD()) {
141 loader->LoadDigits();
142 }
abaf1f1d 143
f7336fa3 144 // Read in the digit arrays
88cb7938 145 return (fDigitsManager->ReadDigits(loader->TreeD()));
f7336fa3 146
147}
148
25ca55ce 149//_____________________________________________________________________________
150Bool_t AliTRDclusterizerV1::ReadDigits(TTree *digitsTree)
151{
152 //
153 // Reads the digits arrays from the input tree
154 //
155
156 // Read in the digit arrays
157 return (fDigitsManager->ReadDigits(digitsTree));
158
159}
160
928e9fae 161//_____________________________________________________________________________
a6dd11e9 162Bool_t AliTRDclusterizerV1::ReadDigits(AliRawReader *rawReader)
928e9fae 163{
164 //
165 // Reads the digits arrays from the ddl file
166 //
167
4ab68796 168 AliTRDrawData raw;
4ab68796 169 fDigitsManager = raw.Raw2Digits(rawReader);
928e9fae 170
171 return kTRUE;
172
173}
174
f7336fa3 175//_____________________________________________________________________________
793ff80c 176Bool_t AliTRDclusterizerV1::MakeClusters()
f7336fa3 177{
178 //
179 // Generates the cluster.
180 //
181
6d50f529 182 Int_t row = 0;
183 Int_t col = 0;
184 Int_t time = 0;
185 Int_t icham = 0;
186 Int_t iplan = 0;
187 Int_t isect = 0;
188 Int_t iPad = 0;
189
190 AliTRDdataArrayI *digitsIn;
a6dd11e9 191 AliTRDdataArrayI *tracksIn;
f7336fa3 192
193 // Get the geometry
6d50f529 194 AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader);
40609f3f 195 if (!geo) {
196 AliWarning("Loading default TRD geometry!");
197 geo = new AliTRDgeometry();
198 }
a6dd11e9 199
6d50f529 200 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
201 if (!calibration) {
4fad09c9 202 AliFatal("No AliTRDcalibDB instance available\n");
3551db50 203 return kFALSE;
204 }
205
6d50f529 206 AliTRDSimParam *simParam = AliTRDSimParam::Instance();
207 if (!simParam) {
208 AliError("No AliTRDSimParam instance available\n");
3becff3c 209 return kFALSE;
210 }
211
6d50f529 212 AliTRDRecParam *recParam = AliTRDRecParam::Instance();
213 if (!recParam) {
214 AliError("No AliTRDRecParam instance available\n");
3551db50 215 return kFALSE;
216 }
217
6d50f529 218 AliTRDCommonParam *commonParam = AliTRDCommonParam::Instance();
219 if (!commonParam) {
220 AliError("Could not get common parameters\n");
3551db50 221 return kFALSE;
222 }
f7336fa3 223
8ec071c8 224 // ADC thresholds
a6dd11e9 225 Float_t ADCthreshold = simParam->GetADCthreshold();
3e1a3ad8 226 // Threshold value for the maximum
a6dd11e9 227 Float_t maxThresh = recParam->GetClusMaxThresh();
3e1a3ad8 228 // Threshold value for the digit signal
a6dd11e9 229 Float_t sigThresh = recParam->GetClusSigThresh();
6d50f529 230
56178ff4 231 // Detector wise calibration object for t0
8ec071c8 232 const AliTRDCalDet *calT0Det = calibration->GetT0Det();
233 // Detector wise calibration object for the gain factors
234 const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet();
56178ff4 235
f7336fa3 236 // Iteration limit for unfolding procedure
8230f242 237 const Float_t kEpsilon = 0.01;
8230f242 238 const Int_t kNclus = 3;
239 const Int_t kNsig = 5;
a6dd11e9 240 const Int_t kNdict = AliTRDdigitsManager::kNDict;
241 const Int_t kNtrack = kNdict * kNclus;
3e1a3ad8 242
a5cadd36 243 Int_t iType = 0;
7ad19338 244 Int_t iUnfold = 0;
a5cadd36 245 Double_t ratioLeft = 1.0;
246 Double_t ratioRight = 1.0;
db30bf0f 247
a6dd11e9 248 Int_t iClusterROC = 0;
249
a5cadd36 250 Double_t padSignal[kNsig];
251 Double_t clusterSignal[kNclus];
252 Double_t clusterPads[kNclus];
f7336fa3 253
6d50f529 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 Int_t nTimeTotal = calibration->GetNumberOfTimeBins();
f7336fa3 261
a6dd11e9 262 Int_t dummy[9] = { 0, 0, 0, 0, 0, 0, 0, 0, 0 };
263
6d50f529 264 AliDebug(1,Form("Number of Time Bins = %d.\n",nTimeTotal));
3becff3c 265
3e1a3ad8 266 // Start clustering in every chamber
6d50f529 267 for (icham = chamBeg; icham < chamEnd; icham++) {
268 for (iplan = planBeg; iplan < planEnd; iplan++) {
269 for (isect = sectBeg; isect < sectEnd; isect++) {
f7336fa3 270
a6dd11e9 271 Int_t idet = geo->GetDetector(iplan,icham,isect);
272
273 // Get the digits
274 digitsIn = fDigitsManager->GetDigits(idet);
275 // This is to take care of switched off super modules
276 if (digitsIn->GetNtime() == 0) {
277 continue;
278 }
279 digitsIn->Expand();
280 AliTRDdataArrayI *tracksTmp = fDigitsManager->GetDictionary(idet,0);
281 tracksTmp->Expand();
282
3becff3c 283 Int_t nRowMax = commonParam->GetRowMax(iplan,icham,isect);
284 Int_t nColMax = commonParam->GetColMax(iplan);
f7336fa3 285
6d50f529 286 AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham);
287
56178ff4 288 // Calibration object with pad wise values for t0
8ec071c8 289 AliTRDCalROC *calT0ROC = calibration->GetT0ROC(idet);
290 // Calibration object with pad wise values for the gain factors
291 AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet);
56178ff4 292 // Calibration value for chamber wise t0
8ec071c8 293 Float_t calT0DetValue = calT0Det->GetValue(idet);
294 // Calibration value for chamber wise gain factor
295 Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet);
56178ff4 296
db30bf0f 297 Int_t nClusters = 0;
298 Int_t nClusters2pad = 0;
299 Int_t nClusters3pad = 0;
300 Int_t nClusters4pad = 0;
301 Int_t nClusters5pad = 0;
302 Int_t nClustersLarge = 0;
3e1a3ad8 303
a6dd11e9 304 // Apply the gain and the tail cancelation via digital filter
6d50f529 305 AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow()
306 ,digitsIn->GetNcol()
8ec071c8 307 ,digitsIn->GetNtime());
308 Transform(digitsIn
309 ,digitsOut
310 ,nRowMax,nColMax,nTimeTotal
311 ,ADCthreshold
312 ,calGainFactorROC
313 ,calGainFactorDetValue);
a6dd11e9 314
315 // Input digits are not needed any more
316 digitsIn->Compress(1,0);
6d50f529 317
3e1a3ad8 318 // Loop through the chamber and find the maxima
319 for ( row = 0; row < nRowMax; row++) {
de4b10e5 320 for ( col = 2; col < nColMax; col++) {
3e1a3ad8 321 for (time = 0; time < nTimeTotal; time++) {
322
11dc3a9e 323 Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time));
3e1a3ad8 324
325 // Look for the maximum
db30bf0f 326 if (signalM >= maxThresh) {
6a242715 327
328 Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time));
329 Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
330
6d50f529 331 if ((TMath::Abs(signalL) <= signalM) &&
21ee7de5 332 (TMath::Abs(signalR) < signalM)) {
6493a90f 333 if ((TMath::Abs(signalL) >= sigThresh) ||
334 (TMath::Abs(signalR) >= sigThresh)) {
a6dd11e9 335 // Maximum found, mark the position by a negative signal
336 digitsOut->SetDataUnchecked(row,col-1,time,-signalM);
337 }
3e1a3ad8 338 }
6a242715 339
3e1a3ad8 340 }
6d50f529 341
342 }
343 }
344 }
a6dd11e9 345 tracksTmp->Compress(1,0);
346
347 // The index to the first cluster of a given ROC
348 Int_t firstClusterROC = -1;
349 // The number of cluster in a given ROC
350 Int_t nClusterROC = 0;
3e1a3ad8 351
352 // Now check the maxima and calculate the cluster position
353 for ( row = 0; row < nRowMax ; row++) {
db30bf0f 354 for (time = 0; time < nTimeTotal; time++) {
355 for ( col = 1; col < nColMax-1; col++) {
3e1a3ad8 356
357 // Maximum found ?
a6dd11e9 358 if (digitsOut->GetDataUnchecked(row,col,time) < 0.0) {
f7336fa3 359
8230f242 360 for (iPad = 0; iPad < kNclus; iPad++) {
3e1a3ad8 361 Int_t iPadCol = col - 1 + iPad;
a6dd11e9 362 clusterSignal[iPad] =
363 TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time));
f7336fa3 364 }
365
db30bf0f 366 // Count the number of pads in the cluster
367 Int_t nPadCount = 0;
a6dd11e9 368 Int_t ii;
369 // Look to the left
370 ii = 0;
6d50f529 371 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) {
db30bf0f 372 nPadCount++;
373 ii++;
374 if (col-ii < 0) break;
375 }
a6dd11e9 376 // Look to the right
db30bf0f 377 ii = 0;
6d50f529 378 while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) {
db30bf0f 379 nPadCount++;
380 ii++;
381 if (col+ii+1 >= nColMax) break;
382 }
db30bf0f 383 nClusters++;
384 switch (nPadCount) {
385 case 2:
386 iType = 0;
387 nClusters2pad++;
388 break;
389 case 3:
390 iType = 1;
391 nClusters3pad++;
392 break;
393 case 4:
394 iType = 2;
395 nClusters4pad++;
396 break;
397 case 5:
398 iType = 3;
399 nClusters5pad++;
400 break;
401 default:
402 iType = 4;
403 nClustersLarge++;
404 break;
405 };
406
6d50f529 407 // Look for 5 pad cluster with minimum in the middle
db30bf0f 408 Bool_t fivePadCluster = kFALSE;
6d50f529 409 if (col < (nColMax - 3)) {
3becff3c 410 if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) {
db30bf0f 411 fivePadCluster = kTRUE;
412 }
6d50f529 413 if ((fivePadCluster) && (col < (nColMax - 5))) {
3becff3c 414 if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) {
db30bf0f 415 fivePadCluster = kFALSE;
416 }
417 }
6d50f529 418 if ((fivePadCluster) && (col > 1)) {
3becff3c 419 if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) {
db30bf0f 420 fivePadCluster = kFALSE;
421 }
422 }
423 }
424
425 // 5 pad cluster
426 // Modify the signal of the overlapping pad for the left part
427 // of the cluster which remains from a previous unfolding
428 if (iUnfold) {
429 clusterSignal[0] *= ratioLeft;
7ad19338 430 iType = 5;
db30bf0f 431 iUnfold = 0;
432 }
433
434 // Unfold the 5 pad cluster
435 if (fivePadCluster) {
436 for (iPad = 0; iPad < kNsig; iPad++) {
3becff3c 437 padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row
a6dd11e9 438 ,col-1+iPad
439 ,time));
f7336fa3 440 }
db30bf0f 441 // Unfold the two maxima and set the signal on
442 // the overlapping pad to the ratio
17b26de4 443 ratioRight = Unfold(kEpsilon,iplan,padSignal);
db30bf0f 444 ratioLeft = 1.0 - ratioRight;
445 clusterSignal[2] *= ratioRight;
7ad19338 446 iType = 5;
db30bf0f 447 iUnfold = 1;
f7336fa3 448 }
f7336fa3 449
a5cadd36 450 Double_t clusterCharge = clusterSignal[0]
451 + clusterSignal[1]
452 + clusterSignal[2];
3e1a3ad8 453
db30bf0f 454 // The position of the cluster
a6dd11e9 455 clusterPads[0] = row + 0.5;
3e1a3ad8 456 // Take the shift of the additional time bins into account
dde59437 457 clusterPads[2] = time + 0.5;
3e1a3ad8 458
3551db50 459 if (recParam->LUTOn()) {
db30bf0f 460 // Calculate the position of the cluster by using the
461 // lookup table method
3becff3c 462 clusterPads[1] = recParam->LUTposition(iplan,clusterSignal[0]
6d50f529 463 ,clusterSignal[1]
464 ,clusterSignal[2]);
db30bf0f 465 }
466 else {
db30bf0f 467 // Calculate the position of the cluster by using the
468 // center of gravity method
a6dd11e9 469 for (Int_t i = 0; i < kNsig; i++) {
470 padSignal[i] = 0.0;
6d50f529 471 }
a6dd11e9 472 padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); // Central pad
473 padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // Left pad
474 padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // Right pad
6d50f529 475 if ((col > 2) &&
476 (TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)) < padSignal[1])) {
3becff3c 477 padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time));
7ad19338 478 }
6d50f529 479 if ((col < nColMax - 3) &&
480 (TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)) < padSignal[3])) {
3becff3c 481 padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time));
7ad19338 482 }
6d50f529 483 clusterPads[1] = GetCOG(padSignal);
db30bf0f 484 }
485
a5cadd36 486 Double_t q0 = clusterSignal[0];
487 Double_t q1 = clusterSignal[1];
488 Double_t q2 = clusterSignal[2];
a6dd11e9 489 Double_t clusterSigmaY2 = (q1 * (q0 + q2) + 4.0 * q0 * q2)
490 / (clusterCharge*clusterCharge);
a819a5f7 491
6d50f529 492 //
3551db50 493 // Calculate the position and the error
6d50f529 494 //
495
56178ff4 496 // Correct for t0 (sum of chamber and pad wise values !!!)
497 Float_t calT0ROCValue = calT0ROC->GetValue(col,row);
498 Int_t clusterTimeBin = TMath::Nint(time - (calT0DetValue + calT0ROCValue));
6d50f529 499 Double_t colSize = padPlane->GetColSize(col);
500 Double_t rowSize = padPlane->GetRowSize(row);
cb86ff6e 501
a5cadd36 502 Double_t clusterPos[3];
a6dd11e9 503 clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1] + 0.5) * colSize;
504 clusterPos[1] = padPlane->GetRowPos(row) - 0.5 * rowSize;
6d50f529 505 clusterPos[2] = CalcXposFromTimebin(clusterPads[2],idet,col,row);
a5cadd36 506 Double_t clusterSig[2];
a6dd11e9 507 clusterSig[0] = (clusterSigmaY2 + 1.0/12.0) * colSize*colSize;
508 clusterSig[1] = rowSize * rowSize / 12.0;
3551db50 509
a6dd11e9 510 // Add the cluster to the output array
511 // The track indices will be stored later
512 AliTRDcluster *cluster = AddCluster(clusterPos
513 ,clusterTimeBin
514 ,idet
515 ,clusterCharge
516 ,dummy
517 ,clusterSig
518 ,iType
524fc8fa 519 ,col
520 ,clusterPads[1]);
521
a6dd11e9 522 // Store the amplitudes of the pads in the cluster for later analysis
523 Short_t signals[7] = { 0, 0, 0, 0, 0, 0, 0 };
6d50f529 524 for (Int_t jPad = col-3; jPad <= col+3; jPad++) {
a6dd11e9 525 if ((jPad < 0) ||
526 (jPad >= nColMax-1)) {
6d50f529 527 continue;
528 }
11dc3a9e 529 signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time)));
c85a4951 530 }
531 cluster->SetSignals(signals);
6d50f529 532
a6dd11e9 533 // Temporarily store the row, column and time bin of the center pad
534 // Used to later on assign the track indices
535 cluster->SetLabel( row,0);
536 cluster->SetLabel( col,1);
537 cluster->SetLabel(time,2);
6d50f529 538
a6dd11e9 539 // Store the index of the first cluster in the current ROC
540 if (firstClusterROC < 0) {
541 firstClusterROC = RecPoints()->GetEntriesFast() - 1;
542 }
543 // Count the number of cluster in the current ROC
544 nClusterROC++;
545
546 } // if: Maximum found ?
547
548 } // loop: pad columns
549 } // loop: time bins
550 } // loop: pad rows
11dc3a9e 551
552 delete digitsOut;
f7336fa3 553
a6dd11e9 554 //
555 // Add the track indices to the found clusters
556 //
557
558 // Temporary array to collect the track indices
559 Int_t *idxTracks = new Int_t[kNtrack*nClusterROC];
560
561 // Loop through the dictionary arrays one-by-one
562 // to keep memory consumption low
563 for (Int_t iDict = 0; iDict < kNdict; iDict++) {
564
565 tracksIn = fDigitsManager->GetDictionary(idet,iDict);
566 tracksIn->Expand();
567
568 // Loop though the clusters found in this ROC
569 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
570
571 AliTRDcluster *cluster = (AliTRDcluster *)
572 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
573 row = cluster->GetLabel(0);
574 col = cluster->GetLabel(1);
575 time = cluster->GetLabel(2);
576
577 for (iPad = 0; iPad < kNclus; iPad++) {
578 Int_t iPadCol = col - 1 + iPad;
579 Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1;
580 idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index;
581 }
582
583 }
584
585 // Compress the arrays
586 tracksIn->Compress(1,0);
587
588 }
589
590 // Copy the track indices into the cluster
591 // Loop though the clusters found in this ROC
592 for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) {
593
594 AliTRDcluster *cluster = (AliTRDcluster *)
595 RecPoints()->UncheckedAt(firstClusterROC+iClusterROC);
596 cluster->SetLabel(-9999,0);
597 cluster->SetLabel(-9999,1);
598 cluster->SetLabel(-9999,2);
599
600 cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]);
601
602 }
603
604 delete [] idxTracks;
f7336fa3 605
3e1a3ad8 606 // Write the cluster and reset the array
793ff80c 607 WriteClusters(idet);
bdbb05bb 608 ResetRecPoints();
6d50f529 609
a6dd11e9 610 } // loop: Sectors
611 } // loop: Planes
612 } // loop: Chambers
f7336fa3 613
f7336fa3 614 return kTRUE;
615
616}
617
a305677e 618//_____________________________________________________________________________
7ad19338 619Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5])
620{
621 //
6d50f529 622 // Get COG position
623 // Used for clusters with more than 3 pads - where LUT not applicable
624 //
625
a6dd11e9 626 Double_t sum = signal[0]
627 + signal[1]
628 + signal[2]
629 + signal[3]
630 + signal[4];
631
632 Double_t res = (0.0 * (-signal[0] + signal[4])
633 + (-signal[1] + signal[3])) / sum;
6d50f529 634
7ad19338 635 return res;
6d50f529 636
7ad19338 637}
638
f7336fa3 639//_____________________________________________________________________________
a6dd11e9 640Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t *padSignal)
f7336fa3 641{
642 //
643 // Method to unfold neighbouring maxima.
644 // The charge ratio on the overlapping pad is calculated
645 // until there is no more change within the range given by eps.
646 // The resulting ratio is then returned to the calling method.
647 //
648
a6dd11e9 649 AliTRDcalibDB *calibration = AliTRDcalibDB::Instance();
6d50f529 650 if (!calibration) {
651 AliError("No AliTRDcalibDB instance available\n");
6a739e92 652 return kFALSE;
3551db50 653 }
6a739e92 654
a5cadd36 655 Int_t irc = 0;
a6dd11e9 656 Int_t itStep = 0; // Count iteration steps
f7336fa3 657
a6dd11e9 658 Double_t ratio = 0.5; // Start value for ratio
659 Double_t prevRatio = 0.0; // Store previous ratio
f7336fa3 660
a6dd11e9 661 Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal
662 Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal
663 Double_t newSignal[3] = { 0.0, 0.0, 0.0 };
f7336fa3 664
3e1a3ad8 665 // Start the iteration
f7336fa3 666 while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) {
667
668 itStep++;
669 prevRatio = ratio;
670
3e1a3ad8 671 // Cluster position according to charge ratio
a5cadd36 672 Double_t maxLeft = (ratio*padSignal[2] - padSignal[0])
673 / (padSignal[0] + padSignal[1] + ratio*padSignal[2]);
674 Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2])
a6dd11e9 675 / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]);
f7336fa3 676
3e1a3ad8 677 // Set cluster charge ratio
6a739e92 678 irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal);
a5cadd36 679 Double_t ampLeft = padSignal[1] / newSignal[1];
6a739e92 680 irc = calibration->PadResponse(1.0,maxRight,plane,newSignal);
a5cadd36 681 Double_t ampRight = padSignal[3] / newSignal[1];
f7336fa3 682
3e1a3ad8 683 // Apply pad response to parameters
6a739e92 684 irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal );
685 irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal);
f7336fa3 686
3e1a3ad8 687 // Calculate new overlapping ratio
a5cadd36 688 ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] /
a6dd11e9 689 (newLeftSignal[2] + newRightSignal[0]));
f7336fa3 690
691 }
692
693 return ratio;
694
695}
696
3becff3c 697//_____________________________________________________________________________
a6dd11e9 698void AliTRDclusterizerV1::Transform(AliTRDdataArrayI *digitsIn
699 , AliTRDdataArrayF *digitsOut
8ec071c8 700 , Int_t nRowMax, Int_t nColMax, Int_t nTimeTotal
701 , Float_t ADCthreshold
702 , AliTRDCalROC *calGainFactorROC
703 , Float_t calGainFactorDetValue)
3becff3c 704{
3becff3c 705 //
cb86ff6e 706 // Apply gain factor
a6dd11e9 707 // Apply tail cancelation: Transform digitsIn to digitsOut
3becff3c 708 //
709
6d50f529 710 Int_t iRow = 0;
711 Int_t iCol = 0;
712 Int_t iTime = 0;
3becff3c 713
a6dd11e9 714 AliTRDRecParam *recParam = AliTRDRecParam::Instance();
6d50f529 715 if (!recParam) {
716 AliError("No AliTRDRecParam instance available\n");
3becff3c 717 return;
718 }
11dc3a9e 719
a6dd11e9 720 Double_t *inADC = new Double_t[nTimeTotal]; // ADC data before tail cancellation
721 Double_t *outADC = new Double_t[nTimeTotal]; // ADC data after tail cancellation
3becff3c 722
6d50f529 723 for (iRow = 0; iRow < nRowMax; iRow++ ) {
724 for (iCol = 0; iCol < nColMax; iCol++ ) {
a6dd11e9 725
56178ff4 726 Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow);
727 Double_t gain = calGainFactorDetValue
728 * calGainFactorROCValue;
729
6d50f529 730 for (iTime = 0; iTime < nTimeTotal; iTime++) {
3becff3c 731
cb86ff6e 732 //
6d50f529 733 // Add gain
cb86ff6e 734 //
a6dd11e9 735 inADC[iTime] = digitsIn->GetDataUnchecked(iRow,iCol,iTime);
736 inADC[iTime] /= gain;
737 outADC[iTime] = inADC[iTime];
6d50f529 738
3becff3c 739 }
740
741 // Apply the tail cancelation via the digital filter
6d50f529 742 if (recParam->TCOn()) {
a305677e 743 DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp());
3becff3c 744 }
745
6d50f529 746 for (iTime = 0; iTime < nTimeTotal; iTime++) {
747
3becff3c 748 // Store the amplitude of the digit if above threshold
a305677e 749 if (outADC[iTime] > ADCthreshold) {
11dc3a9e 750 digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]);
3becff3c 751 }
752
753 }
754
755 }
3becff3c 756 }
757
758 delete [] inADC;
759 delete [] outADC;
760
761 return;
762
763}
764
3becff3c 765//_____________________________________________________________________________
a6dd11e9 766void AliTRDclusterizerV1::DeConvExp(Double_t *source, Double_t *target
767 , Int_t n, Int_t nexp)
3becff3c 768{
769 //
6d50f529 770 // Tail cancellation by deconvolution for PASA v4 TRF
3becff3c 771 //
772
773 Double_t rates[2];
774 Double_t coefficients[2];
775
6d50f529 776 // Initialization (coefficient = alpha, rates = lambda)
3becff3c 777 Double_t R1 = 1.0;
778 Double_t R2 = 1.0;
779 Double_t C1 = 0.5;
780 Double_t C2 = 0.5;
781
782 if (nexp == 1) { // 1 Exponentials
783 R1 = 1.156;
784 R2 = 0.130;
785 C1 = 0.066;
786 C2 = 0.000;
787 }
788 if (nexp == 2) { // 2 Exponentials
789 R1 = 1.156;
790 R2 = 0.130;
791 C1 = 0.114;
792 C2 = 0.624;
793 }
794
795 coefficients[0] = C1;
796 coefficients[1] = C2;
797
a6dd11e9 798 Double_t Dt = 0.1;
3becff3c 799
800 rates[0] = TMath::Exp(-Dt/(R1));
801 rates[1] = TMath::Exp(-Dt/(R2));
802
6d50f529 803 Int_t i = 0;
804 Int_t k = 0;
3becff3c 805
6d50f529 806 Double_t reminder[2];
807 Double_t correction;
808 Double_t result;
3becff3c 809
6d50f529 810 // Attention: computation order is important
811 correction = 0.0;
812 for (k = 0; k < nexp; k++) {
813 reminder[k] = 0.0;
814 }
815 for (i = 0; i < n; i++) {
a6dd11e9 816 result = (source[i] - correction); // No rescaling
3becff3c 817 target[i] = result;
818
6d50f529 819 for (k = 0; k < nexp; k++) {
820 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
821 }
822 correction = 0.0;
823 for (k = 0; k < nexp; k++) {
824 correction += reminder[k];
825 }
3becff3c 826 }
827
828}