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.4 2006/08/28 17:12:18 cblume
19 Remove the last print statements
21 Revision 1.3 2006/08/11 17:58:05 cblume
22 Next round of effc++ changes
24 Revision 1.2 2006/04/05 12:45:40 hristov
25 Updated TRD trigger code. Now the trigger code can run:
27 - in the simulation step, through the central trigger interface, to
28 produce and store "tracklets" and to produce and analyze tracks, with
29 inputs to the central trigger
31 - in the reconstruction step: if the tracklets have already been produced
32 in the simulation step (by the trigger option), then only the tracks are
33 produced and stored in the ESD event; if not, the trigger start by
34 producing the tracklets, etc.
38 Revision 1.1.1.1 2004/08/19 14:58:11 vulpescu
41 Revision 1.1.1.1 2004/08/18 07:47:17 vulpescu
46 ///////////////////////////////////////////////////////////////////////////////
49 // Multi Chip Module exponential filter and tracklet finder //
52 ///////////////////////////////////////////////////////////////////////////////
58 #include "AliTRDmcm.h"
59 #include "AliTRDtrigParam.h"
63 //_____________________________________________________________________________
64 AliTRDmcm::AliTRDmcm()
85 // AliTRDmcm default constructor
91 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
95 for (i = 0; i < kMcmCol; i++) {
97 for (j = 0; j < kMcmTBmax; j++) {
99 fIsClus[i][j] = kFALSE;
105 //_____________________________________________________________________________
106 AliTRDmcm::AliTRDmcm(const Int_t id)
127 // AliTRDmcm constructor
133 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
137 for (i = 0; i < kMcmCol; i++) {
139 for (j = 0; j < kMcmTBmax; j++) {
141 fIsClus[i][j] = kFALSE;
145 fTime1 = AliTRDtrigParam::Instance()->GetTime1();
146 fTime2 = AliTRDtrigParam::Instance()->GetTime2();
147 fClusThr = AliTRDtrigParam::Instance()->GetClusThr();
148 fPadThr = AliTRDtrigParam::Instance()->GetPadThr();
150 AliTRDtrigParam::Instance()->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
154 //_____________________________________________________________________________
155 AliTRDmcm::AliTRDmcm(const AliTRDmcm &m)
161 ,fColFirst(m.fColFirst)
162 ,fColLast(m.fColLast)
165 ,fClusThr(m.fClusThr)
167 ,fNtrkSeeds(m.fNtrkSeeds)
172 ,fPedestal(m.fPedestal)
176 // AliTRDmcm copy constructor
182 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
183 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
184 ((AliTRDmcm &) m).fSeedCol[i] = -1;
186 for (i = 0; i < kMcmCol; i++) {
187 ((AliTRDmcm &) m).fPadHits[i] = 0;
188 for (j = 0; j < kMcmTBmax; j++) {
189 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
190 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
196 //_____________________________________________________________________________
197 AliTRDmcm::~AliTRDmcm()
200 // AliTRDmcm destructor
205 //_____________________________________________________________________________
206 AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
209 // Assignment operator
212 if (this != &m) ((AliTRDmcm &) m).Copy(*this);
218 //_____________________________________________________________________________
219 void AliTRDmcm::Copy(TObject &m) const
228 ((AliTRDmcm &) m).fNtrk = fNtrk;
229 ((AliTRDmcm &) m).fRobId = fRobId;
230 ((AliTRDmcm &) m).fChaId = fChaId;
231 ((AliTRDmcm &) m).fRow = fRow;
232 ((AliTRDmcm &) m).fColFirst = fColFirst;
233 ((AliTRDmcm &) m).fColLast = fColLast;
234 ((AliTRDmcm &) m).fTime1 = fTime1;
235 ((AliTRDmcm &) m).fTime2 = fTime2;
236 ((AliTRDmcm &) m).fClusThr = fClusThr;
237 ((AliTRDmcm &) m).fPadThr = fPadThr;
238 ((AliTRDmcm &) m).fNtrkSeeds = fNtrkSeeds;
239 ((AliTRDmcm &) m).fR1 = fR1;
240 ((AliTRDmcm &) m).fR2 = fR2;
241 ((AliTRDmcm &) m).fC1 = fC1;
242 ((AliTRDmcm &) m).fC2 = fC2;
243 ((AliTRDmcm &) m).fPedestal = fPedestal;
244 ((AliTRDmcm &) m).fId = fId;
246 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
247 ((AliTRDmcm &) m).fTrkIndex[i] = 0;
248 ((AliTRDmcm &) m).fSeedCol[i] = -1;
250 for (i = 0; i < kMcmCol; i++) {
251 ((AliTRDmcm &) m).fPadHits[i] = 0;
252 for (j = 0; j < kMcmTBmax; j++) {
253 ((AliTRDmcm &) m).fADC[i][j] = 0.0;
254 ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
260 //_____________________________________________________________________________
261 void AliTRDmcm::AddTrk(const Int_t id)
264 // Add a tracklet index
267 fTrkIndex[fNtrk] = id;
274 //_____________________________________________________________________________
275 void AliTRDmcm::Reset()
284 for (i = 0; i < kMcmCol; i++) {
286 for (j = 0; j < kMcmTBmax; j++) {
288 fIsClus[i][j] = kFALSE;
291 for (i = 0; i < kMaxTrackletsPerMCM; i++) {
297 //_____________________________________________________________________________
298 Bool_t AliTRDmcm::Run()
304 AliDebug(2,(Form("Run MCM %d\n",Id())));
312 Float_t amp[3] = { 0.0, 0.0, 0.0 };
314 Int_t clusCol[kMcmCol/2];
315 Float_t clusAmp[kMcmCol/2];
320 for (iTime = fTime1; iTime <= fTime2; iTime++) {
324 for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
325 amp[0] = fADC[iCol-1][iTime];
326 amp[1] = fADC[iCol ][iTime];
327 amp[2] = fADC[iCol+1][iTime];
328 if (IsCluster(amp)) {
329 fIsClus[iCol][iTime] = kTRUE;
330 clusCol[nClus] = iCol;
331 clusAmp[nClus] = amp[0]+amp[1]+amp[2];
333 if (nClus == kMcmCol/2) {
334 AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
340 // ...but no more than six...
341 if (nClus > (Int_t) kSelClus) {
342 for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
343 fIsClus[clusCol[j]][iTime] = kFALSE;
347 // ...and take the largest four.
348 Int_t nClusPlus = nClus - kMaxClus;
349 for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
351 for (i = 0; i < nClus; i++) {
352 if (fIsClus[clusCol[i]][iTime]) {
353 if (clusAmp[i] <= veryLarge) {
354 veryLarge = clusAmp[i];
359 fIsClus[clusCol[clusMin]][iTime] = kFALSE;
364 } // end main TB loop
366 if ((fNtrkSeeds = CreateSeeds())) {
374 //_____________________________________________________________________________
375 Int_t AliTRDmcm::CreateSeeds()
378 // Make column seeds (from Falk Lesser, ex KIP)
385 AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
387 // Working array for hit sums
388 Int_t fHit2padSum[2][kMcmCol];
390 // Initialize the array
391 for (i = 0; i < 2; i++) {
392 for (j = 0; j < kMcmCol; j++ ) {
394 fHit2padSum[i][j] = j;
396 fHit2padSum[i][j] = -1;
401 Int_t sum10 = AliTRDtrigParam::Instance()->GetSum10();
402 Int_t sum12 = AliTRDtrigParam::Instance()->GetSum12();
406 for (i = 0; i < kMcmCol; i++) {
407 if (i < (kMcmCol-1)) {
408 if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
409 fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
412 fHit2padSum[1][i] = -1;
416 if (fPadHits[i] >= sum12) {
417 fHit2padSum[1][i] = fPadHits[i];
420 fHit2padSum[1][i] = -1;
423 if (fHit2padSum[1][i] > 0) {
428 // sort the sums in decreasing order of the amplitude
429 Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
431 // arrange (maximum number of) candidates in increasing order of the column number
432 nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
433 Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
435 for (i = 0; i < nSeeds; i++) {
436 fSeedCol[i] = fHit2padSum[0][i];
439 // reject multiple found tracklets
440 Int_t imax = nSeeds - 1;
441 for (i = 0; i < imax; i++) {
443 if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
445 if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
446 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
450 AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
461 //_____________________________________________________________________________
462 void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir)
465 // Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
466 // in the direction: dir = 0 ascending order
467 // dir = 1 descending order
479 for (i = 0; i < (nel-1); i++) {
480 if (x2[i+1] < x2[i]) {
490 } while (sort == kFALSE);
498 for (i = 0; i < (nel-1); i++) {
499 if (x2[i+1] > x2[i]) {
509 } while (sort == kFALSE);
515 //_____________________________________________________________________________
516 void AliTRDmcm::AddTimeBin(const Int_t iTime)
519 // Build column seeds
522 for (Int_t iPad = 1; iPad < (kMcmCol-1); iPad++) {
523 if (fIsClus[iPad][iTime]) {
530 //_____________________________________________________________________________
531 Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
534 // Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
538 if (amp[0] > amp[1] || amp[2] > amp[1]) {
542 // -> cluster amplitude
543 if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
548 if (amp[0] < fPadThr && amp[2] < fPadThr) {
556 //_____________________________________________________________________________
557 void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
560 // Exponential filter
566 Double_t sour[kMcmTBmax];
567 Double_t dtarg[kMcmTBmax];
568 Int_t itarg[kMcmTBmax];
574 for (iCol = 0; iCol < kMcmCol; iCol++) {
575 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
576 sour[iTime] = fADC[iCol][iTime];
578 DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
579 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
580 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
587 for (iCol = 0; iCol < kMcmCol; iCol++) {
588 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
589 sour[iTime] = fADC[iCol][iTime];
591 DeConvExpD(sour,itarg,kMcmTBmax,nexp);
592 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
593 fADC[iCol][iTime] = itarg[iTime];
600 for (iCol = 0; iCol < kMcmCol; iCol++) {
601 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
602 sour[iTime] = fADC[iCol][iTime];
604 DeConvExpMI(sour,dtarg,kMcmTBmax);
605 for (iTime = 0; iTime < kMcmTBmax; iTime++) {
606 fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
613 AliError(Form("Invalid filter type %d ! \n",ftype));
620 //_____________________________________________________________________________
621 void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
624 // Exponential filter "analog"
629 Double_t reminder[2];
633 Double_t coefficients[2];
635 // Initialize (coefficient = alpha, rates = lambda)
637 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
638 Double_t r1, r2, c1, c2;
644 coefficients[0] = c1;
645 coefficients[1] = c2;
648 rates[0] = TMath::Exp(-dt/(r1));
649 rates[1] = TMath::Exp(-dt/(r2));
651 // Attention: computation order is important
653 for (k = 0; k < nexp; k++) {
657 for (i = 0; i < n; i++) {
659 result = (source[i] - correction); // no rescaling
662 for (k = 0; k < nexp; k++) {
663 reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
667 for (k = 0; k < nexp; k++) {
668 correction += reminder[k];
675 //_____________________________________________________________________________
676 void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
679 // Exponential filter "digital"
695 // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
696 // initialize (coefficient = alpha, rates = lambda)
718 fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
719 fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
720 iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
721 iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
724 fAlphaL = (Int_t) (c1 * 2048.0);
725 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
728 fAlphaL = (Int_t) (c1 * 2048.0);
729 fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
730 iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
731 iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
734 Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
735 Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
736 Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
737 Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
745 Int_t iFactor = ((Int_t) fPedestal) << 2;
747 Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
748 rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
749 rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
751 // further initialization
752 if ((rem1 + rem2) > 0x0FFF) {
756 correction = (rem1 + rem2) & 0x0FFF;
759 fTailPed = iFactor - correction;
761 for (i = 0; i < n; i++) {
763 result = ((Int_t)source[i] - correction);
770 h1 = (rem1 + ((iAlphaL * result) >> 11));
778 h2 = (rem2 + ((iAlphaS * result) >> 11));
786 rem1 = (iLambdaL * h1 ) >> 11;
787 rem2 = (iLambdaS * h2 ) >> 11;
789 if ((rem1 + rem2) > 0x0FFF) {
793 correction = (rem1 + rem2) & 0x0FFF;
800 //_____________________________________________________________________________
801 void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
804 // Exponential filter (M. Ivanov)
813 for (i = 0; i < n; i++) {
819 Float_t lambda0 = (1.0 / fR2) * dt;
820 Float_t lambda1 = (1.0 / fR1) * dt;
822 TailMakerSpline(sig1,sig2,lambda0,n);
823 TailCancelationMI(sig2,sig3,0.7,lambda1,n);
825 for (i = 0; i < n; i++) {
831 //______________________________________________________________________________
832 void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
835 // Special filter (M. Ivanov)
840 Double_t l = TMath::Exp(-lambda*0.5);
844 // Initialize in[] and out[] goes 0 ... 2*n+19
845 for (i = 0; i < n*2+20; i++) {
852 in[1] = (ampin[0] + ampin[1]) * 0.5;
854 // Add charge to the end
855 for (i = 0; i < 22; i++) {
856 // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
857 in[2*(n-1)+i] = ampin[n-1];
860 // Use arithmetic mean
861 for (i = 1; i < n-1; i++) {
862 // in[] goes 2, 3, ... , 2*n-4, 2*n-3
864 in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
870 for (i = 2*n; i >= 0; i--) {
871 out[i] = in[i] + temp;
872 temp = l*(temp+in[i]);
875 for (i = 0; i < n; i++){
876 //ampout[i] = out[2*i+1]; // org
877 ampout[i] = out[2*i];
882 //______________________________________________________________________________
883 void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
884 , Double_t norm, Double_t lambda, Int_t n)
887 // Special filter (M. Ivanov)
892 Double_t l = TMath::Exp(-lambda*0.5);
893 Double_t k = l*(1.0 - norm*lambda*0.5);
897 // Initialize in[] and out[] goes 0 ... 2*n+19
898 for (i = 0; i < n*2+20; i++) {
905 in[1] = (ampin[0]+ampin[1])*0.5;
907 // Add charge to the end
908 for (i =-2; i < 22; i++) {
909 // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
910 in[2*(n-1)+i] = ampin[n-1];
913 for (i = 1; i < n-2; i++) {
914 // in[] goes 2, 3, ... , 2*n-6, 2*n-5
916 in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
917 //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
923 for (i = 1; i <= 2*n; i++) {
924 out[i] = in[i] + (k-l)*temp;
925 temp = in[i] + k *temp;
928 for (i = 0; i < n; i++) {
929 //ampout[i] = out[2*i+1]; // org
930 //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
931 ampout[i] = TMath::Max(out[2*i],0.0);