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.3 2006/08/11 17:58:05 cblume
19 Next round of effc++ changes
21 Revision 1.2 2006/04/05 12:45:40 hristov
22 Updated TRD trigger code. Now the trigger code can run:
24 - in the simulation step, through the central trigger interface, to
25 produce and store "tracklets" and to produce and analyze tracks, with
26 inputs to the central trigger
28 - in the reconstruction step: if the tracklets have already been produced
29 in the simulation step (by the trigger option), then only the tracks are
30 produced and stored in the ESD event; if not, the trigger start by
31 producing the tracklets, etc.
35 Revision 1.1.1.1 2004/08/19 14:58:11 vulpescu
38 Revision 1.1.1.1 2004/08/18 07:47:17 vulpescu
43 ///////////////////////////////////////////////////////////////////////////////
46 // Multi Chip Module exponential filter and tracklet finder //
49 ///////////////////////////////////////////////////////////////////////////////
55 #include "AliTRDmcm.h"
56 #include "AliTRDtrigParam.h"
60 //_____________________________________________________________________________
61 AliTRDmcm::AliTRDmcm()
83 // AliTRDmcm default constructor
89 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
93 for (i = 0; i < kMcmCol; i++) {
95 for (j = 0; j < kMcmTBmax; j++) {
97 fIsClus[i][j] = kFALSE;
103 //_____________________________________________________________________________
104 AliTRDmcm::AliTRDmcm(AliTRDtrigParam *trigp, const Int_t id)
113 ,fTime1(trigp->GetTime1())
114 ,fTime2(trigp->GetTime2())
115 ,fClusThr(trigp->GetClusThr())
116 ,fPadThr(trigp->GetPadThr())
126 // AliTRDmcm constructor
132 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
136 for (i = 0; i < kMcmCol; i++) {
138 for (j = 0; j < kMcmTBmax; j++) {
140 fIsClus[i][j] = kFALSE;
144 fTrigParam->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
148 //_____________________________________________________________________________
149 AliTRDmcm::AliTRDmcm(const AliTRDmcm &m)
156 ,fColFirst(m.fColFirst)
157 ,fColLast(m.fColLast)
160 ,fClusThr(m.fClusThr)
162 ,fNtrkSeeds(m.fNtrkSeeds)
167 ,fPedestal(m.fPedestal)
171 // AliTRDmcm copy constructor
177 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
178 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
179 ((AliTRDmcm &) m).fSeedCol[i] = -1;
181 for (i = 0; i < kMcmCol; i++) {
182 ((AliTRDmcm &) m).fPadHits[i] = 0;
183 for (j = 0; j < kMcmTBmax; j++) {
184 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
185 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
191 //_____________________________________________________________________________
192 AliTRDmcm::~AliTRDmcm()
195 // AliTRDmcm destructor
200 //_____________________________________________________________________________
201 AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
204 // Assignment operator
207 if (this != &m) ((AliTRDmcm &) m).Copy(*this);
213 //_____________________________________________________________________________
214 void AliTRDmcm::Copy(TObject &m) const
223 ((AliTRDmcm &) m).fTrigParam = NULL;
224 ((AliTRDmcm &) m).fNtrk = fNtrk;
225 ((AliTRDmcm &) m).fRobId = fRobId;
226 ((AliTRDmcm &) m).fChaId = fChaId;
227 ((AliTRDmcm &) m).fRow = fRow;
228 ((AliTRDmcm &) m).fColFirst = fColFirst;
229 ((AliTRDmcm &) m).fColLast = fColLast;
230 ((AliTRDmcm &) m).fTime1 = fTime1;
231 ((AliTRDmcm &) m).fTime2 = fTime2;
232 ((AliTRDmcm &) m).fClusThr = fClusThr;
233 ((AliTRDmcm &) m).fPadThr = fPadThr;
234 ((AliTRDmcm &) m).fNtrkSeeds = fNtrkSeeds;
235 ((AliTRDmcm &) m).fR1 = fR1;
236 ((AliTRDmcm &) m).fR2 = fR2;
237 ((AliTRDmcm &) m).fC1 = fC1;
238 ((AliTRDmcm &) m).fC2 = fC2;
239 ((AliTRDmcm &) m).fPedestal = fPedestal;
240 ((AliTRDmcm &) m).fId = fId;
242 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
243 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
244 ((AliTRDmcm &) m).fSeedCol[i] = -1;
246 for (i = 0; i < kMcmCol; i++) {
247 ((AliTRDmcm &) m).fPadHits[i] = 0;
248 for (j = 0; j < kMcmTBmax; j++) {
249 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
250 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
256 //_____________________________________________________________________________
257 void AliTRDmcm::AddTrk(const Int_t id)
260 // Add a tracklet index
263 fTrkIndex[fNtrk] = id;
270 //_____________________________________________________________________________
271 void AliTRDmcm::Reset()
280 for (i = 0; i < kMcmCol; i++) {
282 for (j = 0; j < kMcmTBmax; j++) {
284 fIsClus[i][j] = kFALSE;
287 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
293 //_____________________________________________________________________________
294 Bool_t AliTRDmcm::Run()
300 AliDebug(2,(Form("Run MCM %d\n",Id())));
309 Float_t amp[3] = { 0.0, 0.0, 0.0 };
311 Int_t clusCol[kMcmCol/2];
312 Float_t clusAmp[kMcmCol/2];
317 for (iTime = fTime1; iTime <= fTime2; iTime++) {
321 for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
322 amp[0] = fADC[iCol-1][iTime];
323 amp[1] = fADC[iCol ][iTime];
324 amp[2] = fADC[iCol+1][iTime];
325 if (IsCluster(amp)) {
326 fIsClus[iCol][iTime] = kTRUE;
327 clusCol[nClus] = iCol;
328 clusAmp[nClus] = amp[0]+amp[1]+amp[2];
330 if (nClus == kMcmCol/2) {
331 AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
337 // ...but no more than six...
338 if (nClus > (Int_t) kSelClus) {
339 for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
340 fIsClus[clusCol[j]][iTime] = kFALSE;
344 // ...and take the largest four.
345 Int_t nClusPlus = nClus - kMaxClus;
346 for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
348 for (i = 0; i < nClus; i++) {
349 if (fIsClus[clusCol[i]][iTime]) {
350 if (clusAmp[i] <= veryLarge) {
351 veryLarge = clusAmp[i];
356 fIsClus[clusCol[clusMin]][iTime] = kFALSE;
361 } // end main TB loop
363 if ((fNtrkSeeds = CreateSeeds())) {
371 //_____________________________________________________________________________
372 Int_t AliTRDmcm::CreateSeeds()
375 // Make column seeds (from Falk Lesser, ex KIP)
382 AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
384 // Working array for hit sums
385 Int_t fHit2padSum[2][kMcmCol];
387 // Initialize the array
388 for (i = 0; i < 2; i++) {
389 for (j = 0; j < kMcmCol; j++ ) {
391 fHit2padSum[i][j] = j;
393 fHit2padSum[i][j] = -1;
398 Int_t sum10 = fTrigParam->GetSum10();
399 Int_t sum12 = fTrigParam->GetSum12();
403 for (i = 0; i < kMcmCol; i++) {
404 if (i < (kMcmCol-1)) {
405 if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
406 fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
409 fHit2padSum[1][i] = -1;
413 if (fPadHits[i] >= sum12) {
414 fHit2padSum[1][i] = fPadHits[i];
417 fHit2padSum[1][i] = -1;
420 if (fHit2padSum[1][i] > 0) {
425 // sort the sums in decreasing order of the amplitude
426 Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
428 // arrange (maximum number of) candidates in increasing order of the column number
429 nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
430 Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
432 for (i = 0; i < nSeeds; i++) {
433 fSeedCol[i] = fHit2padSum[0][i];
436 // reject multiple found tracklets
437 Int_t imax = nSeeds - 1;
438 for (i = 0; i < imax; i++) {
440 if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
442 if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
443 if (fTrigParam->GetDebugLevel() > 1) {
444 AliWarning(Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
449 if (fTrigParam->GetDebugLevel() > 1) {
450 AliWarning(Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
462 //_____________________________________________________________________________
463 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir)
466 // Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
467 // in the direction: dir = 0 ascending order
468 // dir = 1 descending order
480 for (i = 0; i < (nel-1); i++) {
481 if (x2[i+1] < x2[i]) {
491 } while (sort == kFALSE);
499 for (i = 0; i < (nel-1); i++) {
500 if (x2[i+1] > x2[i]) {
510 } while (sort == kFALSE);
516 //_____________________________________________________________________________
517 void AliTRDmcm::AddTimeBin(const Int_t iTime)
520 // Build column seeds
523 for (Int_t iPad = 1; iPad < (kMcmCol-1); iPad++) {
524 if (fIsClus[iPad][iTime]) {
531 //_____________________________________________________________________________
532 Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
535 // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
539 if (amp[0] > amp[1] || amp[2] > amp[1]) {
543 // -> cluster amplitude
544 if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
549 if (amp[0] < fPadThr && amp[2] < fPadThr) {
557 //_____________________________________________________________________________
558 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
561 // Exponential filter
567 Double_t sour[kMcmTBmax];
568 Double_t dtarg[kMcmTBmax];
569 Int_t itarg[kMcmTBmax];
575 for (iCol = 0; iCol < kMcmCol; iCol++) {
576 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
577 sour[iTime] = fADC[iCol][iTime];
579 DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
580 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
581 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
588 for (iCol = 0; iCol < kMcmCol; iCol++) {
589 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
590 sour[iTime] = fADC[iCol][iTime];
592 DeConvExpD(sour,itarg,kMcmTBmax,nexp);
593 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
594 fADC[iCol][iTime] = itarg[iTime];
601 for (iCol = 0; iCol < kMcmCol; iCol++) {
602 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
603 sour[iTime] = fADC[iCol][iTime];
605 DeConvExpMI(sour,dtarg,kMcmTBmax);
606 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
607 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
614 AliError(Form("Invalid filter type %d ! \n",ftype));
621 //_____________________________________________________________________________
622 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
625 // Exponential filter "analog"
630 Double_t reminder[2];
634 Double_t coefficients[2];
636 // Initialize (coefficient = alpha, rates = lambda)
638 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
639 Double_t r1, r2, c1, c2;
645 coefficients[0] = c1;
646 coefficients[1] = c2;
649 rates[0] = TMath::Exp(-dt/(r1));
650 rates[1] = TMath::Exp(-dt/(r2));
652 // Attention: computation order is important
654 for (k = 0; k < nexp; k++) {
658 for (i = 0; i < n; i++) {
660 result = (source[i] - correction); // no rescaling
663 for (k = 0; k < nexp; k++) {
664 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
668 for (k = 0; k < nexp; k++) {
669 correction += reminder[k];
676 //_____________________________________________________________________________
677 void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
680 // Exponential filter "digital"
696 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
697 // initialize (coefficient = alpha, rates = lambda)
719 fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
720 fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
721 iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
722 iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
725 fAlphaL = (Int_t) (c1 * 2048.0);
726 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
729 fAlphaL = (Int_t) (c1 * 2048.0);
730 fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
731 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
732 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
735 Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
736 Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
737 Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
738 Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
746 Int_t iFactor = ((Int_t) fPedestal) << 2;
748 Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
749 rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
750 rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
752 // further initialization
753 if ((rem1 + rem2) > 0x0FFF) {
757 correction = (rem1 + rem2) & 0x0FFF;
760 fTailPed = iFactor - correction;
762 for (i = 0; i < n; i++) {
764 result = ((Int_t)source[i] - correction);
771 h1 = (rem1 + ((iAlphaL * result) >> 11));
779 h2 = (rem2 + ((iAlphaS * result) >> 11));
787 rem1 = (iLambdaL * h1 ) >> 11;
788 rem2 = (iLambdaS * h2 ) >> 11;
790 if ((rem1 + rem2) > 0x0FFF) {
794 correction = (rem1 + rem2) & 0x0FFF;
801 //_____________________________________________________________________________
802 void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
805 // Exponential filter (M. Ivanov)
814 for (i = 0; i < n; i++) {
820 Float_t lambda0 = (1.0 / fR2) * dt;
821 Float_t lambda1 = (1.0 / fR1) * dt;
823 TailMakerSpline(sig1,sig2,lambda0,n);
824 TailCancelationMI(sig2,sig3,0.7,lambda1,n);
826 for (i = 0; i < n; i++) {
832 //______________________________________________________________________________
833 void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
836 // Special filter (M. Ivanov)
841 Double_t l = TMath::Exp(-lambda*0.5);
845 // Initialize in[] and out[] goes 0 ... 2*n+19
846 for (i = 0; i < n*2+20; i++) {
853 in[1] = (ampin[0] + ampin[1]) * 0.5;
855 // Add charge to the end
856 for (i = 0; i < 22; i++) {
857 // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
858 in[2*(n-1)+i] = ampin[n-1];
861 // Use arithmetic mean
862 for (i = 1; i < n-1; i++) {
863 // in[] goes 2, 3, ... , 2*n-4, 2*n-3
865 in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
871 for (i = 2*n; i >= 0; i--) {
872 out[i] = in[i] + temp;
873 temp = l*(temp+in[i]);
876 for (i = 0; i < n; i++){
877 //ampout[i] = out[2*i+1]; // org
878 ampout[i] = out[2*i];
883 //______________________________________________________________________________
884 void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
885 , Double_t norm, Double_t lambda, Int_t n)
888 // Special filter (M. Ivanov)
893 Double_t l = TMath::Exp(-lambda*0.5);
894 Double_t k = l*(1.0 - norm*lambda*0.5);
898 // Initialize in[] and out[] goes 0 ... 2*n+19
899 for (i = 0; i < n*2+20; i++) {
906 in[1] = (ampin[0]+ampin[1])*0.5;
908 // Add charge to the end
909 for (i =-2; i < 22; i++) {
910 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
911 in[2*(n-1)+i] = ampin[n-1];
914 for (i = 1; i < n-2; i++) {
915 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
917 in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
918 //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
924 for (i = 1; i <= 2*n; i++) {
925 out[i] = in[i] + (k-l)*temp;
926 temp = in[i] + k *temp;
929 for (i = 0; i < n; i++) {
930 //ampout[i] = out[2*i+1]; // org
931 //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
932 ampout[i] = TMath::Max(out[2*i],0.0);