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.1.1.1 2004/08/19 14:58:11 vulpescu
21 Revision 1.1.1.1 2004/08/18 07:47:17 vulpescu
26 ///////////////////////////////////////////////////////////////////////////////
29 // Multi Chip Module exponential filter and tracklet finder //
32 ///////////////////////////////////////////////////////////////////////////////
36 #include "AliTRDmcm.h"
37 #include "AliTRDtrigParam.h"
41 //_____________________________________________________________________________
42 AliTRDmcm::AliTRDmcm()
45 // AliTRDmcm default constructor
51 for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
59 for (Int_t i = 0; i < kMcmCol; i++) {
61 for (Int_t j = 0; j < kMcmTBmax; j++) {
63 fIsClus[i][j] = kFALSE;
71 for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
85 //_____________________________________________________________________________
86 AliTRDmcm::AliTRDmcm(AliTRDtrigParam *trigp, const Int_t id)
89 // AliTRDmcm constructor
95 for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
103 for (Int_t i = 0; i < kMcmCol; i++) {
105 for (Int_t j = 0; j < kMcmTBmax; j++) {
107 fIsClus[i][j] = kFALSE;
110 fPadThr = fTrigParam->GetPadThr();
111 fClusThr = fTrigParam->GetClusThr();
112 fTime1 = fTrigParam->GetTime1();
113 fTime2 = fTrigParam->GetTime2();
115 for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
125 fTrigParam->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
131 //_____________________________________________________________________________
132 AliTRDmcm::~AliTRDmcm()
135 // AliTRDmcm destructor
140 //_____________________________________________________________________________
141 AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
144 // assignment operator
147 if (this != &m) ((AliTRDmcm &) m).Copy(*this);
152 //_____________________________________________________________________________
153 void AliTRDmcm::Copy(TObject &m) const
159 ((AliTRDmcm &) m).fTrigParam = NULL;
160 ((AliTRDmcm &) m).fNtrk = fNtrk;
161 ((AliTRDmcm &) m).fRobId = fRobId;
162 ((AliTRDmcm &) m).fChaId = fChaId;
163 ((AliTRDmcm &) m).fRow = fRow;
164 ((AliTRDmcm &) m).fColFirst = fColFirst;
165 ((AliTRDmcm &) m).fColLast = fColLast;
166 ((AliTRDmcm &) m).fTime1 = fTime1;
167 ((AliTRDmcm &) m).fTime2 = fTime2;
168 ((AliTRDmcm &) m).fClusThr = fClusThr;
169 ((AliTRDmcm &) m).fPadThr = fPadThr;
170 ((AliTRDmcm &) m).fNtrkSeeds = fNtrkSeeds;
171 ((AliTRDmcm &) m).fR1 = fR1;
172 ((AliTRDmcm &) m).fR2 = fR2;
173 ((AliTRDmcm &) m).fC1 = fC1;
174 ((AliTRDmcm &) m).fC2 = fC2;
175 ((AliTRDmcm &) m).fPedestal = fPedestal;
176 ((AliTRDmcm &) m).fId = fId;
178 for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
179 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
181 for (Int_t i = 0; i < kMcmCol; i++) {
182 ((AliTRDmcm &) m).fPadHits[i] = 0;
183 for (Int_t j = 0; j < kMcmTBmax; j++) {
184 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
185 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
188 for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
189 ((AliTRDmcm &) m).fSeedCol[i] = -1;
194 //_____________________________________________________________________________
195 void AliTRDmcm::AddTrk(const Int_t id)
198 // Add a tracklet index
201 fTrkIndex[fNtrk] = id;
208 //_____________________________________________________________________________
209 void AliTRDmcm::Reset()
215 for (Int_t i = 0; i < kMcmCol; i++) {
217 for (Int_t j = 0; j < kMcmTBmax; j++) {
219 fIsClus[i][j] = kFALSE;
222 for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
228 //_____________________________________________________________________________
229 Bool_t AliTRDmcm::Run()
235 if ( fTrigParam->GetDebugLevel() > 1 ) printf("AliTRDmcm::Run MCM %d\n",Id());
237 Float_t amp[3] = {0.0, 0.0, 0.0};
239 Int_t clusCol[kMcmCol/2];
240 Float_t clusAmp[kMcmCol/2];
244 for (Int_t iTime = fTime1; iTime <= fTime2; iTime++) { // main TB loop
248 for (Int_t iCol = 1; iCol < (kMcmCol-1); iCol++) {
249 amp[0] = fADC[iCol-1][iTime];
250 amp[1] = fADC[iCol ][iTime];
251 amp[2] = fADC[iCol+1][iTime];
252 if (IsCluster(amp)) {
253 fIsClus[iCol][iTime] = kTRUE;
254 clusCol[nClus] = iCol;
255 clusAmp[nClus] = amp[0]+amp[1]+amp[2];
257 if (nClus == kMcmCol/2) {
258 printf("Too many clusters in time bin %2d MCM %d...\n",iTime,Id());
265 // ...but no more than six...
266 if (nClus > (Int_t)kSelClus) {
267 for (Int_t j = kSelClus/2; j < nClus-kSelClus/2; j++) {
268 fIsClus[clusCol[j]][iTime] = kFALSE;
272 // ...and take the largest four.
274 Int_t nClusPlus = nClus - kMaxClus;
275 for (Int_t iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
277 for (Int_t i = 0; i < nClus; i++) {
278 if (fIsClus[clusCol[i]][iTime]) {
279 if (clusAmp[i] <= veryLarge) {
280 veryLarge = clusAmp[i];
285 fIsClus[clusCol[clusMin]][iTime] = kFALSE;
290 } // end main TB loop
292 if (fTrigParam->GetDebugLevel() > 1) {
293 for (Int_t i = fTime1; i <= fTime2; i++) {
295 for (Int_t j = 0; j < kMcmCol; j++) {
296 printf("%1d ",fIsClus[j][i]);
301 for (Int_t iPad = 0; iPad < kMcmCol; iPad++) {
302 printf("%2d ",fPadHits[iPad]);
307 if ((fNtrkSeeds = CreateSeeds())) {
316 //_____________________________________________________________________________
317 Int_t AliTRDmcm::CreateSeeds()
320 // Make column seeds (from Falk Lesser, ex KIP)
323 if ( fTrigParam->GetDebugLevel() > 1 ) printf("AliTRDmcm::CreateSeeds MCM %d \n",Id());
327 // working array for hit sums
328 Int_t fHit2padSum[2][kMcmCol];
330 // initialize the array
331 for( Int_t i = 0; i < 2; i++ ) {
332 for( Int_t j = 0; j < kMcmCol; j++ ) {
334 fHit2padSum[i][j] = j;
336 fHit2padSum[i][j] = -1;
341 Int_t sum10 = fTrigParam->GetSum10();
342 Int_t sum12 = fTrigParam->GetSum12();
346 for( Int_t i = 0; i < kMcmCol; i++ ) {
347 if( i < (kMcmCol-1) ) {
348 if( (fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12) ) {
349 fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
351 fHit2padSum[1][i] = -1;
354 if ( fPadHits[i] >= sum12 ) {
355 fHit2padSum[1][i] = fPadHits[i];
357 fHit2padSum[1][i] = -1;
360 if (fHit2padSum[1][i] > 0) nsum2seed++;
363 if (fTrigParam->GetDebugLevel() > 1) {
364 printf("fHit2padSum: ");
365 for( Int_t i = 0; i < kMcmCol; i++ ) {
366 printf("%2d ",fHit2padSum[0][i]);
370 for( Int_t i = 0; i < kMcmCol; i++ ) {
371 printf("%2d ",fHit2padSum[1][i]);
376 // sort the sums in decreasing order of the amplitude
377 Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
379 if (fTrigParam->GetDebugLevel() > 1) {
380 printf("fHit2padSum: ");
381 for( Int_t i = 0; i < kMcmCol; i++ ) {
382 printf("%2d ",fHit2padSum[0][i]);
386 for( Int_t i = 0; i < kMcmCol; i++ ) {
387 printf("%2d ",fHit2padSum[1][i]);
392 // arrange (maximum number of) candidates in increasing order of the column number
393 nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
394 Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
396 for (Int_t i = 0; i < nSeeds; i++) {
397 fSeedCol[i] = fHit2padSum[0][i];
400 if (fTrigParam->GetDebugLevel() > 1) {
401 printf("Found %d seeds before multiple rejection. \n",nSeeds);
402 printf("fHit2padSum: ");
403 for( Int_t i = 0; i < kMcmCol; i++ ) {
404 printf("%2d ",fHit2padSum[0][i]);
408 for( Int_t i = 0; i < kMcmCol; i++ ) {
409 printf("%2d ",fHit2padSum[1][i]);
414 // reject multiple found tracklets
415 Int_t imax = nSeeds-1;
416 for (Int_t i = 0; i < imax; i++) {
418 if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
420 if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
421 if (fTrigParam->GetDebugLevel() > 1)
422 printf("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]);
425 if (fTrigParam->GetDebugLevel() > 1)
426 printf("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]);
433 if ( fTrigParam->GetDebugLevel() > 1 ) {
434 printf("Found %d seeds in MCM %d ",nSeeds,Id());
435 for (Int_t i = 0; i < (imax+1); i++) {
436 if (fSeedCol[i] >= 0) printf(", %02d ",fSeedCol[i]);
445 //_____________________________________________________________________________
446 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir)
449 // Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
450 // in the direction: dir = 0 ascending order
451 // dir = 1 descending order
461 for ( Int_t i = 0; i < (nel-1); i++ )
462 if ( x2[i+1] < x2[i] ) {
471 } while ( sort == kFALSE );
479 for ( Int_t i = 0; i < (nel-1); i++ )
480 if ( x2[i+1] > x2[i] ) {
489 } while ( sort == kFALSE );
495 //_____________________________________________________________________________
496 void AliTRDmcm::AddTimeBin(const Int_t iTime)
499 // Build column seeds
502 for (Int_t iPad = 1; iPad < (kMcmCol-1); iPad++) {
503 if (fIsClus[iPad][iTime]) {
510 //_____________________________________________________________________________
511 Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
514 // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
518 if (amp[0] > amp[1] || amp[2] > amp[1]) return kFALSE;
520 // -> cluster amplitude
521 if ((amp[0]+amp[1]+amp[2]) < fClusThr) return kFALSE;
524 if (amp[0] < fPadThr && amp[2] < fPadThr) return kFALSE;
530 //_____________________________________________________________________________
531 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
534 // exponential filter
537 Double_t sour[kMcmTBmax];
538 Double_t dtarg[kMcmTBmax];
539 Int_t itarg[kMcmTBmax];
545 for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
546 for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
547 sour[iTime] = fADC[iCol][iTime];
549 DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
550 for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
551 //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,dtarg[iTime]);
552 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
559 for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
560 for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
561 sour[iTime] = fADC[iCol][iTime];
563 DeConvExpD(sour,itarg,kMcmTBmax,nexp);
564 for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
565 fADC[iCol][iTime] = itarg[iTime];
572 for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
573 for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
574 sour[iTime] = fADC[iCol][iTime];
576 DeConvExpMI(sour,dtarg,kMcmTBmax);
577 for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
578 //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,dtarg[iTime]);
579 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
586 printf("Invalid filter type %d ! \n",ftype);
593 //_____________________________________________________________________________
594 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
597 // Exponential filter "analog"
601 Double_t coefficients[2];
603 // initialize (coefficient = alpha, rates = lambda)
605 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
606 Double_t r1, r2, c1, c2;
612 coefficients[0] = c1;
613 coefficients[1] = c2;
616 rates[0] = TMath::Exp(-dt/(r1));
617 rates[1] = TMath::Exp(-dt/(r2));
620 Double_t reminder[2];
621 Double_t correction, result;
623 /* attention: computation order is important */
625 for ( k=0; k<nexp; k++ ) reminder[k]=0.0;
627 for ( i=0; i<n; i++ ) {
628 result = ( source[i] - correction ); // no rescaling
631 for ( k=0; k<nexp; k++ ) reminder[k] = rates[k] * ( reminder[k] + coefficients[k] * result);
634 for ( k=0; k<nexp; k++ ) correction += reminder[k];
639 //_____________________________________________________________________________
640 void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
643 // Exponential filter "digital"
646 Int_t fAlphaL, fAlphaS, fLambdaL, fLambdaS, fTailPed;
647 Int_t iAlphaL, iAlphaS, iLambdaL, iLambdaS;
649 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
650 // initialize (coefficient = alpha, rates = lambda)
652 fLambdaL = 0; fAlphaL = 0; fLambdaS = 0; fAlphaS = 0;
653 iLambdaL = 0; iAlphaL = 0; iLambdaS = 0; iAlphaS = 0;
657 Double_t r1, r2, c1, c2;
663 fLambdaL = (Int_t)((TMath::Exp(-dt/r1)-0.75)*2048.0);
664 fLambdaS = (Int_t)((TMath::Exp(-dt/r2)-0.25)*2048.0);
665 iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
666 iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
669 fAlphaL = (Int_t)(c1*2048.0);
670 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
673 fAlphaL = (Int_t)(c1*2048.0);
674 fAlphaS = (Int_t)((c2-0.5)*2048.0);
675 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
676 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
679 Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
680 Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
681 Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
682 Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
686 Int_t correction, result;
687 Int_t iFactor = ((Int_t)fPedestal)<<2;
689 Double_t xi = 1-(iLl*iAs + iLs*iAl); // calculation of equilibrium values of the
690 rem1 = (Int_t)((iFactor/xi) * ((1-iLs)*iLl*iAl)); // internal registers to prevent switch on effects.
691 rem2 = (Int_t)((iFactor/xi) * ((1-iLl)*iLs*iAs));
693 // further initialization
694 if ( (rem1 + rem2) > 0x0FFF ) {
697 correction = (rem1 + rem2) & 0x0FFF;
700 fTailPed = iFactor - correction;
702 for (Int_t i=0; i<n; i++) {
704 result = ( (Int_t)source[i] - correction );
711 h1 = ( rem1 + ((iAlphaL * result) >> 11) );
718 h2 = ( rem2 + ((iAlphaS * result) >> 11));
725 rem1 = (iLambdaL * h1 ) >> 11;
726 rem2 = (iLambdaS * h2 ) >> 11;
728 if ( (rem1 + rem2) > 0x0FFF ) {
731 correction = (rem1 + rem2) & 0x0FFF;
737 //_____________________________________________________________________________
738 void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
741 // Exponential filter (M. Ivanov)
744 Double_t sig1[100], sig2[100], sig3[100];//, sig4[100];
745 for (Int_t i = 0; i < n; i++) sig1[i] = source[i];
749 //Float_t lambda0 = 9.8016*dt; // short
750 //Float_t lambda1 = 1.0778*dt; // long
752 Float_t lambda0 = (1.0/fR2)*dt;
753 Float_t lambda1 = (1.0/fR1)*dt;
755 TailMakerSpline(sig1,sig2,lambda0,n);
756 TailCancelationMI(sig2,sig3,0.7,lambda1,n);
758 for (Int_t i = 0; i < n; i++) target[i] = sig3[i];
762 //______________________________________________________________________________
763 void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
766 // Special filter (M. Ivanov)
769 Double_t l = TMath::Exp(-lambda*0.5);
774 // initialize in[] and out[] goes 0 ... 2*n+19
775 for (Int_t i=0; i<n*2+20; i++){
782 in[1] = (ampin[0]+ampin[1])*0.5;
784 // add charge to the end
785 for (Int_t i=0; i<22; i++){
786 // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
787 in[2*(n-1)+i] = ampin[n-1];
791 for (Int_t i=1; i<n-1; i++){
792 // in[] goes 2, 3, ... , 2*n-4, 2*n-3
794 in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
797 // add charge to the end
798 for (Int_t i=-2; i<22; i++){
799 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
800 in[2*(n-1)+i] = ampin[n-1];
804 for (Int_t i=1; i<n-2; i++){
805 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
807 in[2*i+1] = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;
814 for (int i=2*n; i>=0; i--){
815 out[i] = in[i] + temp;
816 temp = l*(temp+in[i]);
820 for (int i=0;i<n;i++){
821 //ampout[i] = out[2*i+1]; // org
822 ampout[i] = out[2*i];
827 //______________________________________________________________________________
828 void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout, Double_t norm, Double_t lambda, Int_t n)
831 // Special filter (M. Ivanov)
834 Double_t l = TMath::Exp(-lambda*0.5);
835 Double_t k = l*(1.-norm*lambda*0.5);
840 // initialize in[] and out[] goes 0 ... 2*n+19
841 for (Int_t i=0; i<n*2+20; i++){
848 in[1] = (ampin[0]+ampin[1])*0.5;
850 // add charge to the end
851 for (Int_t i=-2; i<22; i++){
852 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
853 in[2*(n-1)+i] = ampin[n-1];
857 for (Int_t i=1; i<n-2; i++){
858 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
860 in[2*i+1] = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;
861 //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
868 for (int i=1; i<=2*n; i++){
869 out[i] = in[i] + (k-l)*temp;
870 temp = in[i] + k*temp;
874 for (int i=0; i<n; i++){
875 //ampout[i] = out[2*i+1]; // org
876 //ampout[i] = TMath::Max(out[2*i+1], 0.0); // org
877 ampout[i] = TMath::Max(out[2*i], 0.0);