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 ////////////////////////////////////////////////////////////////////////////
21 // Multi Chip Module exponential filter and tracklet finder //
24 ////////////////////////////////////////////////////////////////////////////
30 #include "AliTRDmcm.h"
31 #include "AliTRDtrigParam.h"
35 //_____________________________________________________________________________
36 AliTRDmcm::AliTRDmcm()
57 // AliTRDmcm default constructor
63 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
67 for (i = 0; i < kMcmCol; i++) {
69 for (j = 0; j < kMcmTBmax; j++) {
71 fIsClus[i][j] = kFALSE;
77 //_____________________________________________________________________________
78 AliTRDmcm::AliTRDmcm(const Int_t id)
99 // AliTRDmcm constructor
105 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
109 for (i = 0; i < kMcmCol; i++) {
111 for (j = 0; j < kMcmTBmax; j++) {
113 fIsClus[i][j] = kFALSE;
117 fTime1 = AliTRDtrigParam::Instance()->GetTime1();
118 fTime2 = AliTRDtrigParam::Instance()->GetTime2();
119 fClusThr = AliTRDtrigParam::Instance()->GetClusThr();
120 fPadThr = AliTRDtrigParam::Instance()->GetPadThr();
122 AliTRDtrigParam::Instance()->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
126 //_____________________________________________________________________________
127 AliTRDmcm::AliTRDmcm(const AliTRDmcm &m)
133 ,fColFirst(m.fColFirst)
134 ,fColLast(m.fColLast)
137 ,fClusThr(m.fClusThr)
139 ,fNtrkSeeds(m.fNtrkSeeds)
144 ,fPedestal(m.fPedestal)
148 // AliTRDmcm copy constructor
154 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
155 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
156 ((AliTRDmcm &) m).fSeedCol[i] = -1;
158 for (i = 0; i < kMcmCol; i++) {
159 ((AliTRDmcm &) m).fPadHits[i] = 0;
160 for (j = 0; j < kMcmTBmax; j++) {
161 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
162 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
168 //_____________________________________________________________________________
169 AliTRDmcm::~AliTRDmcm()
172 // AliTRDmcm destructor
177 //_____________________________________________________________________________
178 AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
181 // Assignment operator
184 if (this != &m) ((AliTRDmcm &) m).Copy(*this);
190 //_____________________________________________________________________________
191 void AliTRDmcm::Copy(TObject &m) const
200 ((AliTRDmcm &) m).fNtrk = fNtrk;
201 ((AliTRDmcm &) m).fRobId = fRobId;
202 ((AliTRDmcm &) m).fChaId = fChaId;
203 ((AliTRDmcm &) m).fRow = fRow;
204 ((AliTRDmcm &) m).fColFirst = fColFirst;
205 ((AliTRDmcm &) m).fColLast = fColLast;
206 ((AliTRDmcm &) m).fTime1 = fTime1;
207 ((AliTRDmcm &) m).fTime2 = fTime2;
208 ((AliTRDmcm &) m).fClusThr = fClusThr;
209 ((AliTRDmcm &) m).fPadThr = fPadThr;
210 ((AliTRDmcm &) m).fNtrkSeeds = fNtrkSeeds;
211 ((AliTRDmcm &) m).fR1 = fR1;
212 ((AliTRDmcm &) m).fR2 = fR2;
213 ((AliTRDmcm &) m).fC1 = fC1;
214 ((AliTRDmcm &) m).fC2 = fC2;
215 ((AliTRDmcm &) m).fPedestal = fPedestal;
216 ((AliTRDmcm &) m).fId = fId;
218 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
219 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
220 ((AliTRDmcm &) m).fSeedCol[i] = -1;
222 for (i = 0; i < kMcmCol; i++) {
223 ((AliTRDmcm &) m).fPadHits[i] = 0;
224 for (j = 0; j < kMcmTBmax; j++) {
225 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
226 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
232 //_____________________________________________________________________________
233 void AliTRDmcm::AddTrk(const Int_t id)
236 // Add a tracklet index
239 fTrkIndex[fNtrk] = id;
246 //_____________________________________________________________________________
247 void AliTRDmcm::Reset()
256 for (i = 0; i < kMcmCol; i++) {
258 for (j = 0; j < kMcmTBmax; j++) {
260 fIsClus[i][j] = kFALSE;
263 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
269 //_____________________________________________________________________________
270 Bool_t AliTRDmcm::Run()
276 AliDebug(2,(Form("Run MCM %d\n",Id())));
284 Float_t amp[3] = { 0.0, 0.0, 0.0 };
286 Int_t clusCol[kMcmCol/2];
287 Float_t clusAmp[kMcmCol/2];
292 for (iTime = fTime1; iTime <= fTime2; iTime++) {
296 for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
297 amp[0] = fADC[iCol-1][iTime];
298 amp[1] = fADC[iCol ][iTime];
299 amp[2] = fADC[iCol+1][iTime];
300 if (IsCluster(amp)) {
301 fIsClus[iCol][iTime] = kTRUE;
302 clusCol[nClus] = iCol;
303 clusAmp[nClus] = amp[0]+amp[1]+amp[2];
305 if (nClus == kMcmCol/2) {
306 AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
312 // ...but no more than six...
313 if (nClus > (Int_t) kSelClus) {
314 for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
315 fIsClus[clusCol[j]][iTime] = kFALSE;
319 // ...and take the largest four.
320 Int_t nClusPlus = nClus - kMaxClus;
321 for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
323 for (i = 0; i < nClus; i++) {
324 if (fIsClus[clusCol[i]][iTime]) {
325 if (clusAmp[i] <= veryLarge) {
326 veryLarge = clusAmp[i];
331 fIsClus[clusCol[clusMin]][iTime] = kFALSE;
336 } // end main TB loop
338 if ((fNtrkSeeds = CreateSeeds())) {
346 //_____________________________________________________________________________
347 Int_t AliTRDmcm::CreateSeeds()
350 // Make column seeds (from Falk Lesser, ex KIP)
357 AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
359 // Working array for hit sums
360 Int_t fHit2padSum[2][kMcmCol];
362 // Initialize the array
363 for (i = 0; i < 2; i++) {
364 for (j = 0; j < kMcmCol; j++ ) {
366 fHit2padSum[i][j] = j;
368 fHit2padSum[i][j] = -1;
373 Int_t sum10 = AliTRDtrigParam::Instance()->GetSum10();
374 Int_t sum12 = AliTRDtrigParam::Instance()->GetSum12();
378 for (i = 0; i < kMcmCol; i++) {
379 if (i < (kMcmCol-1)) {
380 if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
381 fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
384 fHit2padSum[1][i] = -1;
388 if (fPadHits[i] >= sum12) {
389 fHit2padSum[1][i] = fPadHits[i];
392 fHit2padSum[1][i] = -1;
395 if (fHit2padSum[1][i] > 0) {
400 // sort the sums in decreasing order of the amplitude
401 Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
403 // arrange (maximum number of) candidates in increasing order of the column number
404 nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
405 Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
407 for (i = 0; i < nSeeds; i++) {
408 fSeedCol[i] = fHit2padSum[0][i];
411 // reject multiple found tracklets
412 Int_t imax = nSeeds - 1;
413 for (i = 0; i < imax; i++) {
415 if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
417 if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
418 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
422 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
433 //_____________________________________________________________________________
434 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir) const
437 // Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
438 // in the direction: dir = 0 ascending order
439 // dir = 1 descending order
451 for (i = 0; i < (nel-1); i++) {
452 if (x2[i+1] < x2[i]) {
462 } while (sort == kFALSE);
470 for (i = 0; i < (nel-1); i++) {
471 if (x2[i+1] > x2[i]) {
481 } while (sort == kFALSE);
487 //_____________________________________________________________________________
488 void AliTRDmcm::AddTimeBin(const Int_t iTime)
491 // Build column seeds
494 for (Int_t iPad = 1; iPad < (kMcmCol-1); iPad++) {
495 if (fIsClus[iPad][iTime]) {
502 //_____________________________________________________________________________
503 Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
506 // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
510 if (amp[0] > amp[1] || amp[2] > amp[1]) {
514 // -> cluster amplitude
515 if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
520 if (amp[0] < fPadThr && amp[2] < fPadThr) {
528 //_____________________________________________________________________________
529 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
532 // Exponential filter
538 Double_t sour[kMcmTBmax];
539 Double_t dtarg[kMcmTBmax];
540 Int_t itarg[kMcmTBmax];
546 for (iCol = 0; iCol < kMcmCol; iCol++) {
547 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
548 sour[iTime] = fADC[iCol][iTime];
550 DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
551 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
552 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
559 for (iCol = 0; iCol < kMcmCol; iCol++) {
560 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
561 sour[iTime] = fADC[iCol][iTime];
563 DeConvExpD(sour,itarg,kMcmTBmax,nexp);
564 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
565 fADC[iCol][iTime] = itarg[iTime];
572 for (iCol = 0; iCol < kMcmCol; iCol++) {
573 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
574 sour[iTime] = fADC[iCol][iTime];
576 DeConvExpMI(sour,dtarg,kMcmTBmax);
577 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
578 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
585 AliError(Form("Invalid filter type %d ! \n",ftype));
592 //_____________________________________________________________________________
593 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
596 // Exponential filter "analog"
601 Double_t reminder[2];
605 Double_t coefficients[2];
607 // Initialize (coefficient = alpha, rates = lambda)
609 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
610 Double_t r1, r2, c1, c2;
616 coefficients[0] = c1;
617 coefficients[1] = c2;
620 rates[0] = TMath::Exp(-dt/(r1));
621 rates[1] = TMath::Exp(-dt/(r2));
623 // Attention: computation order is important
625 for (k = 0; k < nexp; k++) {
629 for (i = 0; i < n; i++) {
631 result = (source[i] - correction); // no rescaling
634 for (k = 0; k < nexp; k++) {
635 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
639 for (k = 0; k < nexp; k++) {
640 correction += reminder[k];
647 //_____________________________________________________________________________
648 void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
651 // Exponential filter "digital"
667 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
668 // initialize (coefficient = alpha, rates = lambda)
690 fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
691 fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
692 iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
693 iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
696 fAlphaL = (Int_t) (c1 * 2048.0);
697 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
700 fAlphaL = (Int_t) (c1 * 2048.0);
701 fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
702 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
703 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
706 Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
707 Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
708 Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
709 Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
717 Int_t iFactor = ((Int_t) fPedestal) << 2;
719 Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
720 rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
721 rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
723 // further initialization
724 if ((rem1 + rem2) > 0x0FFF) {
728 correction = (rem1 + rem2) & 0x0FFF;
731 fTailPed = iFactor - correction;
733 for (i = 0; i < n; i++) {
735 result = ((Int_t)source[i] - correction);
742 h1 = (rem1 + ((iAlphaL * result) >> 11));
750 h2 = (rem2 + ((iAlphaS * result) >> 11));
758 rem1 = (iLambdaL * h1 ) >> 11;
759 rem2 = (iLambdaS * h2 ) >> 11;
761 if ((rem1 + rem2) > 0x0FFF) {
765 correction = (rem1 + rem2) & 0x0FFF;
772 //_____________________________________________________________________________
773 void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
776 // Exponential filter (M. Ivanov)
785 for (i = 0; i < n; i++) {
791 Float_t lambda0 = (1.0 / fR2) * dt;
792 Float_t lambda1 = (1.0 / fR1) * dt;
794 TailMakerSpline(sig1,sig2,lambda0,n);
795 TailCancelationMI(sig2,sig3,0.7,lambda1,n);
797 for (i = 0; i < n; i++) {
803 //______________________________________________________________________________
804 void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
807 // Special filter (M. Ivanov)
812 Double_t l = TMath::Exp(-lambda*0.5);
816 // Initialize in[] and out[] goes 0 ... 2*n+19
817 for (i = 0; i < n*2+20; i++) {
824 in[1] = (ampin[0] + ampin[1]) * 0.5;
826 // Add charge to the end
827 for (i = 0; i < 22; i++) {
828 // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
829 in[2*(n-1)+i] = ampin[n-1];
832 // Use arithmetic mean
833 for (i = 1; i < n-1; i++) {
834 // in[] goes 2, 3, ... , 2*n-4, 2*n-3
836 in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
842 for (i = 2*n; i >= 0; i--) {
843 out[i] = in[i] + temp;
844 temp = l*(temp+in[i]);
847 for (i = 0; i < n; i++){
848 //ampout[i] = out[2*i+1]; // org
849 ampout[i] = out[2*i];
854 //______________________________________________________________________________
855 void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
856 , Double_t norm, Double_t lambda, Int_t n)
859 // Special filter (M. Ivanov)
864 Double_t l = TMath::Exp(-lambda*0.5);
865 Double_t k = l*(1.0 - norm*lambda*0.5);
869 // Initialize in[] and out[] goes 0 ... 2*n+19
870 for (i = 0; i < n*2+20; i++) {
877 in[1] = (ampin[0]+ampin[1])*0.5;
879 // Add charge to the end
880 for (i =-2; i < 22; i++) {
881 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
882 in[2*(n-1)+i] = ampin[n-1];
885 for (i = 1; i < n-2; i++) {
886 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
888 in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
889 //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
895 for (i = 1; i <= 2*n; i++) {
896 out[i] = in[i] + (k-l)*temp;
897 temp = in[i] + k *temp;
900 for (i = 0; i < n; i++) {
901 //ampout[i] = out[2*i+1]; // org
902 //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
903 ampout[i] = TMath::Max(out[2*i],0.0);