* provided "as is" without express or implied warranty. *
**************************************************************************/
-/*
-$Log$
-Revision 1.1.1.1 2004/08/19 14:58:11 vulpescu
-CVS head
+/* $Id$ */
-Revision 1.1.1.1 2004/08/18 07:47:17 vulpescu
-test
-
-*/
+////////////////////////////////////////////////////////////////////////////
+// //
+// //
+// Multi Chip Module exponential filter and tracklet finder //
+// //
+// //
+////////////////////////////////////////////////////////////////////////////
#include <TMath.h>
-#include <TH2F.h>
+
+#include "AliLog.h"
#include "AliTRDmcm.h"
#include "AliTRDtrigParam.h"
//_____________________________________________________________________________
AliTRDmcm::AliTRDmcm()
+ :TObject()
+ ,fNtrk(0)
+ ,fRobId(0)
+ ,fChaId(0)
+ ,fRow(0)
+ ,fColFirst(0)
+ ,fColLast(0)
+ ,fTime1(0)
+ ,fTime2(0)
+ ,fClusThr(0)
+ ,fPadThr(0)
+ ,fNtrkSeeds(0)
+ ,fR1(0)
+ ,fR2(0)
+ ,fC1(0)
+ ,fC2(0)
+ ,fPedestal(0)
+ ,fId(0)
{
-
//
// AliTRDmcm default constructor
//
- fTrigParam = 0;
+ Int_t i = 0;
+ Int_t j = 0;
- fNtrk = 0;
- for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
+ for (i = 0; i < kMaxTrackletsPerMCM; i++) {
fTrkIndex[i] = 0;
+ fSeedCol[i] = -1;
}
- fRobId = 0;
- fChaId = 0;
- fRow = 0;
- fColFirst = 0;
- fColLast = 0;
- for (Int_t i = 0; i < kMcmCol; i++) {
- fPadHits[i] = 0;
- for (Int_t j = 0; j < kMcmTBmax; j++) {
- fADC[i][j] = 0.0;
+ for (i = 0; i < kMcmCol; i++) {
+ fPadHits[i] = 0;
+ for (j = 0; j < kMcmTBmax; j++) {
+ fADC[i][j] = 0.0;
fIsClus[i][j] = kFALSE;
}
}
- fPadThr = 0;
- fClusThr = 0;
- fTime1 = 0;
- fTime2 = 0;
- fNtrkSeeds = 0;
- for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
- fSeedCol[i] = -1;
- }
-
- fR1 = 0.0;
- fR2 = 0.0;
- fC1 = 0.0;
- fC2 = 0.0;
- fPedestal = 0.0;
-
- fId = 0;
}
//_____________________________________________________________________________
-AliTRDmcm::AliTRDmcm(AliTRDtrigParam *trigp, const Int_t id)
+AliTRDmcm::AliTRDmcm(Int_t id)
+ :TObject()
+ ,fNtrk(0)
+ ,fRobId(0)
+ ,fChaId(0)
+ ,fRow(0)
+ ,fColFirst(0)
+ ,fColLast(0)
+ ,fTime1(0)
+ ,fTime2(0)
+ ,fClusThr(0)
+ ,fPadThr(0)
+ ,fNtrkSeeds(0)
+ ,fR1(0)
+ ,fR2(0)
+ ,fC1(0)
+ ,fC2(0)
+ ,fPedestal(0)
+ ,fId(id)
{
//
// AliTRDmcm constructor
//
- fTrigParam = trigp;
+ Int_t i = 0;
+ Int_t j = 0;
- fNtrk = 0;
- for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
+ for (i = 0; i < kMaxTrackletsPerMCM; i++) {
fTrkIndex[i] = 0;
+ fSeedCol[i] = -1;
}
- fRobId = 0;
- fChaId = 0;
- fRow = 0;
- fColFirst = 0;
- fColLast = 0;
- for (Int_t i = 0; i < kMcmCol; i++) {
- fPadHits[i] = 0;
- for (Int_t j = 0; j < kMcmTBmax; j++) {
- fADC[i][j] = 0.0;
+ for (i = 0; i < kMcmCol; i++) {
+ fPadHits[i] = 0;
+ for (j = 0; j < kMcmTBmax; j++) {
+ fADC[i][j] = 0.0;
fIsClus[i][j] = kFALSE;
}
}
- fPadThr = fTrigParam->GetPadThr();
- fClusThr = fTrigParam->GetClusThr();
- fTime1 = fTrigParam->GetTime1();
- fTime2 = fTrigParam->GetTime2();
- fNtrkSeeds = 0;
- for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
- fSeedCol[i] = -1;
- }
+
+ fTime1 = AliTRDtrigParam::Instance()->GetTime1();
+ fTime2 = AliTRDtrigParam::Instance()->GetTime2();
+ fClusThr = AliTRDtrigParam::Instance()->GetClusThr();
+ fPadThr = AliTRDtrigParam::Instance()->GetPadThr();
- fR1 = 0.0;
- fR2 = 0.0;
- fC1 = 0.0;
- fC2 = 0.0;
- fPedestal = 0.0;
+ AliTRDtrigParam::Instance()->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
+
+}
+
+//_____________________________________________________________________________
+AliTRDmcm::AliTRDmcm(const AliTRDmcm &m)
+ :TObject(m)
+ ,fNtrk(m.fNtrk)
+ ,fRobId(m.fRobId)
+ ,fChaId(m.fChaId)
+ ,fRow(m.fRow)
+ ,fColFirst(m.fColFirst)
+ ,fColLast(m.fColLast)
+ ,fTime1(m.fTime1)
+ ,fTime2(m.fTime2)
+ ,fClusThr(m.fClusThr)
+ ,fPadThr(m.fPadThr)
+ ,fNtrkSeeds(m.fNtrkSeeds)
+ ,fR1(m.fR1)
+ ,fR2(m.fR2)
+ ,fC1(m.fC1)
+ ,fC2(m.fC2)
+ ,fPedestal(m.fPedestal)
+ ,fId(m.fId)
+{
+ //
+ // AliTRDmcm copy constructor
+ //
- fTrigParam->GetFilterParam(fR1,fR2,fC1,fC2,fPedestal);
+ Int_t i = 0;
+ Int_t j = 0;
- fId = id;
+ for (i = 0; i < kMaxTrackletsPerMCM; i++) {
+ ((AliTRDmcm &) m).fTrkIndex[i] = 0;
+ ((AliTRDmcm &) m).fSeedCol[i] = -1;
+ }
+ for (i = 0; i < kMcmCol; i++) {
+ ((AliTRDmcm &) m).fPadHits[i] = 0;
+ for (j = 0; j < kMcmTBmax; j++) {
+ ((AliTRDmcm &) m).fADC[i][j] = 0.0;
+ ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
+ }
+ }
}
//_____________________________________________________________________________
AliTRDmcm::~AliTRDmcm()
{
-
//
// AliTRDmcm destructor
//
}
//_____________________________________________________________________________
-void AliTRDmcm::AddTrk(const Int_t id)
+AliTRDmcm &AliTRDmcm::operator=(const AliTRDmcm &m)
+{
+ //
+ // Assignment operator
+ //
+
+ if (this != &m) ((AliTRDmcm &) m).Copy(*this);
+
+ return *this;
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmcm::Copy(TObject &m) const
+{
+ //
+ // Copy function
+ //
+
+ Int_t i = 0;
+ Int_t j = 0;
+
+ ((AliTRDmcm &) m).fNtrk = fNtrk;
+ ((AliTRDmcm &) m).fRobId = fRobId;
+ ((AliTRDmcm &) m).fChaId = fChaId;
+ ((AliTRDmcm &) m).fRow = fRow;
+ ((AliTRDmcm &) m).fColFirst = fColFirst;
+ ((AliTRDmcm &) m).fColLast = fColLast;
+ ((AliTRDmcm &) m).fTime1 = fTime1;
+ ((AliTRDmcm &) m).fTime2 = fTime2;
+ ((AliTRDmcm &) m).fClusThr = fClusThr;
+ ((AliTRDmcm &) m).fPadThr = fPadThr;
+ ((AliTRDmcm &) m).fNtrkSeeds = fNtrkSeeds;
+ ((AliTRDmcm &) m).fR1 = fR1;
+ ((AliTRDmcm &) m).fR2 = fR2;
+ ((AliTRDmcm &) m).fC1 = fC1;
+ ((AliTRDmcm &) m).fC2 = fC2;
+ ((AliTRDmcm &) m).fPedestal = fPedestal;
+ ((AliTRDmcm &) m).fId = fId;
+
+ for (i = 0; i < kMaxTrackletsPerMCM; i++) {
+ ((AliTRDmcm &) m).fTrkIndex[i] = 0;
+ ((AliTRDmcm &) m).fSeedCol[i] = -1;
+ }
+ for (i = 0; i < kMcmCol; i++) {
+ ((AliTRDmcm &) m).fPadHits[i] = 0;
+ for (j = 0; j < kMcmTBmax; j++) {
+ ((AliTRDmcm &) m).fADC[i][j] = 0.0;
+ ((AliTRDmcm &) m).fIsClus[i][j] = kFALSE;
+ }
+ }
+
+}
+
+//_____________________________________________________________________________
+void AliTRDmcm::AddTrk(Int_t id)
{
//
// Add a tracklet index
// Reset MCM data
//
- for (Int_t i = 0; i < kMcmCol; i++) {
+ Int_t i = 0;
+ Int_t j = 0;
+
+ for (i = 0; i < kMcmCol; i++) {
fPadHits[i] = 0;
- for (Int_t j = 0; j < kMcmTBmax; j++) {
- fADC[i][j] = 0.0;
+ for (j = 0; j < kMcmTBmax; j++) {
+ fADC[i][j] = 0.0;
fIsClus[i][j] = kFALSE;
}
}
- for (Int_t i = 0; i < kMaxTrackletsPerMCM; i++) {
+ for (i = 0; i < kMaxTrackletsPerMCM; i++) {
fSeedCol[i] = -1;
}
// Run MCM
//
- if ( fTrigParam->GetDebugLevel() > 1 ) printf("AliTRDmcm::Run MCM %d\n",Id());
+ AliDebug(2,(Form("Run MCM %d\n",Id())));
+
+ Int_t iTime = 0;
+ Int_t iCol = 0;
+ Int_t iPlus = 0;
+ Int_t i = 0;
+ Int_t j = 0;
- Float_t Amp[3] = {0.0, 0.0, 0.0};
+ Float_t amp[3] = { 0.0, 0.0, 0.0 };
Int_t nClus;
- Int_t ClusCol[kMcmCol/2];
- Float_t ClusAmp[kMcmCol/2];
- Float_t VeryLarge;
- Int_t ClusMin = -1;
+ Int_t clusCol[kMcmCol/2];
+ Float_t clusAmp[kMcmCol/2];
+ Float_t veryLarge;
+ Int_t clusMin = -1;
- for (Int_t iTime = fTime1; iTime <= fTime2; iTime++) { // main TB loop
+ // Main TB loop
+ for (iTime = fTime1; iTime <= fTime2; iTime++) {
- // find clusters...
+ // Find clusters...
nClus = 0;
- for (Int_t iCol = 1; iCol < (kMcmCol-1); iCol++) {
- Amp[0] = fADC[iCol-1][iTime];
- Amp[1] = fADC[iCol ][iTime];
- Amp[2] = fADC[iCol+1][iTime];
- if (IsCluster(Amp)) {
+ for (iCol = 1; iCol < (kMcmCol-1); iCol++) {
+ amp[0] = fADC[iCol-1][iTime];
+ amp[1] = fADC[iCol ][iTime];
+ amp[2] = fADC[iCol+1][iTime];
+ if (IsCluster(amp)) {
fIsClus[iCol][iTime] = kTRUE;
- ClusCol[nClus] = iCol;
- ClusAmp[nClus] = Amp[0]+Amp[1]+Amp[2];
+ clusCol[nClus] = iCol;
+ clusAmp[nClus] = amp[0]+amp[1]+amp[2];
nClus++;
if (nClus == kMcmCol/2) {
- printf("Too many clusters in time bin %2d MCM %d...\n",iTime,Id());
- //return kFALSE;
+ AliWarning(Form("Too many clusters in time bin %2d MCM %d...\n",iTime,Id()));
break;
}
}
}
// ...but no more than six...
- if (nClus > (Int_t)kSelClus) {
- for (Int_t j = kSelClus/2; j < nClus-kSelClus/2; j++) {
- fIsClus[ClusCol[j]][iTime] = kFALSE;
+ if (nClus > (Int_t) kSelClus) {
+ for (j = kSelClus/2; j < nClus-kSelClus/2; j++) {
+ fIsClus[clusCol[j]][iTime] = kFALSE;
}
}
// ...and take the largest four.
-
Int_t nClusPlus = nClus - kMaxClus;
- for (Int_t iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
- VeryLarge = 1.E+10;
- for (Int_t i = 0; i < nClus; i++) {
- if (fIsClus[ClusCol[i]][iTime]) {
- if (ClusAmp[i] <= VeryLarge) {
- VeryLarge = ClusAmp[i];
- ClusMin = i;
+ for (iPlus = 0; iPlus < nClusPlus; iPlus++ ) {
+ veryLarge = 1.E+10;
+ for (i = 0; i < nClus; i++) {
+ if (fIsClus[clusCol[i]][iTime]) {
+ if (clusAmp[i] <= veryLarge) {
+ veryLarge = clusAmp[i];
+ clusMin = i;
}
}
}
- fIsClus[ClusCol[ClusMin]][iTime] = kFALSE;
+ fIsClus[clusCol[clusMin]][iTime] = kFALSE;
}
AddTimeBin(iTime);
} // end main TB loop
-
- if (fTrigParam->GetDebugLevel() > 1) {
- for (Int_t i = fTime1; i <= fTime2; i++) {
- printf("%2d: ",i);
- for (Int_t j = 0; j < kMcmCol; j++) {
- printf("%1d ",fIsClus[j][i]);
- }
- printf("\n");
- }
- printf("PadHits: ");
- for (Int_t iPad = 0; iPad < kMcmCol; iPad++) {
- printf("%2d ",fPadHits[iPad]);
- }
- printf("\n");
- }
-
+
if ((fNtrkSeeds = CreateSeeds())) {
-
return kTRUE;
-
}
return kFALSE;
}
+
//_____________________________________________________________________________
Int_t AliTRDmcm::CreateSeeds()
{
// Make column seeds (from Falk Lesser, ex KIP)
//
- if ( fTrigParam->GetDebugLevel() > 1 ) printf("AliTRDmcm::CreateSeeds MCM %d \n",Id());
-
+ Int_t i = 0;
+ Int_t j = 0;
Int_t nSeeds = 0;
- // working array for hit sums
+ AliDebug(2,Form("AliTRDmcm::CreateSeeds MCM %d \n",Id()));
+
+ // Working array for hit sums
Int_t fHit2padSum[2][kMcmCol];
- // initialize the array
- for( Int_t i = 0; i < 2; i++ ) {
- for( Int_t j = 0; j < kMcmCol; j++ ) {
- if( i == 0 ) {
+ // Initialize the array
+ for (i = 0; i < 2; i++) {
+ for (j = 0; j < kMcmCol; j++ ) {
+ if (i == 0) {
fHit2padSum[i][j] = j;
} else {
fHit2padSum[i][j] = -1;
}
}
- Int_t Sum10 = fTrigParam->GetSum10();
- Int_t Sum12 = fTrigParam->GetSum12();
+ Int_t sum10 = AliTRDtrigParam::Instance()->GetSum10();
+ Int_t sum12 = AliTRDtrigParam::Instance()->GetSum12();
- // build the 2padSum
- Int_t Nsum2seed = 0;
- for( Int_t i = 0; i < kMcmCol; i++ ) {
- if( i < (kMcmCol-1) ) {
- if( (fPadHits[i] >= Sum10) && ((fPadHits[i] + fPadHits[i+1]) >= Sum12) ) {
+ // Build the 2padSum
+ Int_t nsum2seed = 0;
+ for (i = 0; i < kMcmCol; i++) {
+ if (i < (kMcmCol-1)) {
+ if ((fPadHits[i] >= sum10) && ((fPadHits[i] + fPadHits[i+1]) >= sum12)) {
fHit2padSum[1][i] = fPadHits[i] + fPadHits[i+1];
- } else {
+ }
+ else {
fHit2padSum[1][i] = -1;
}
- } else {
- if ( fPadHits[i] >= Sum12 ) {
+ }
+ else {
+ if (fPadHits[i] >= sum12) {
fHit2padSum[1][i] = fPadHits[i];
- } else {
+ }
+ else {
fHit2padSum[1][i] = -1;
}
}
- if (fHit2padSum[1][i] > 0) Nsum2seed++;
- }
-
- if (fTrigParam->GetDebugLevel() > 1) {
- printf("fHit2padSum: ");
- for( Int_t i = 0; i < kMcmCol; i++ ) {
- printf("%2d ",fHit2padSum[0][i]);
- }
- printf("\n");
- printf(" ");
- for( Int_t i = 0; i < kMcmCol; i++ ) {
- printf("%2d ",fHit2padSum[1][i]);
+ if (fHit2padSum[1][i] > 0) {
+ nsum2seed++;
}
- printf("\n");
}
// sort the sums in decreasing order of the amplitude
Sort(kMcmCol,&fHit2padSum[0][0],&fHit2padSum[1][0],1);
- if (fTrigParam->GetDebugLevel() > 1) {
- printf("fHit2padSum: ");
- for( Int_t i = 0; i < kMcmCol; i++ ) {
- printf("%2d ",fHit2padSum[0][i]);
- }
- printf("\n");
- printf(" ");
- for( Int_t i = 0; i < kMcmCol; i++ ) {
- printf("%2d ",fHit2padSum[1][i]);
- }
- printf("\n");
- }
-
// arrange (maximum number of) candidates in increasing order of the column number
- nSeeds = TMath::Min(Nsum2seed,kMaxTrackletsPerMCM);
+ nSeeds = TMath::Min(nsum2seed,kMaxTrackletsPerMCM);
Sort(nSeeds,&fHit2padSum[1][0],&fHit2padSum[0][0],0);
- for (Int_t i = 0; i < nSeeds; i++) {
+ for (i = 0; i < nSeeds; i++) {
fSeedCol[i] = fHit2padSum[0][i];
}
- if (fTrigParam->GetDebugLevel() > 1) {
- printf("Found %d seeds before multiple rejection. \n",nSeeds);
- printf("fHit2padSum: ");
- for( Int_t i = 0; i < kMcmCol; i++ ) {
- printf("%2d ",fHit2padSum[0][i]);
- }
- printf("\n");
- printf(" ");
- for( Int_t i = 0; i < kMcmCol; i++ ) {
- printf("%2d ",fHit2padSum[1][i]);
- }
- printf("\n");
- }
-
// reject multiple found tracklets
- Int_t imax = nSeeds-1;
- for (Int_t i = 0; i < imax; i++) {
+ Int_t imax = nSeeds - 1;
+ for (i = 0; i < imax; i++) {
if ((fHit2padSum[0][i]+1) == fHit2padSum[0][i+1]) {
nSeeds--;
if (fHit2padSum[1][i] >= fHit2padSum[1][i+1]) {
- if (fTrigParam->GetDebugLevel() > 1)
- printf("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]);
+ AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i+1]));
fSeedCol[i+1] = -1;
- } else {
- if (fTrigParam->GetDebugLevel() > 1)
- printf("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]);
- fSeedCol[i] = -1;
+ }
+ else {
+ AliDebug(2,Form("Reject seed %1d in col %02d. \n",i,fHit2padSum[0][i]));
+ fSeedCol[i] = -1;
}
}
}
- if ( fTrigParam->GetDebugLevel() > 1 ) {
- printf("Found %d seeds in MCM %d ",nSeeds,Id());
- for (Int_t i = 0; i < (imax+1); i++) {
- if (fSeedCol[i] >= 0) printf(", %02d ",fSeedCol[i]);
- }
- printf("\n");
- }
-
return nSeeds;
}
//_____________________________________________________________________________
-void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir)
+void AliTRDmcm::Sort(Int_t nel, Int_t *x1, Int_t *x2, Int_t dir) const
{
-
+ //
// Sort two parallel vectors (x1[nel], x2[nel]) after the second one (x2)
// in the direction: dir = 0 ascending order
// dir = 1 descending order
+ //
+ Int_t i = 0;
Bool_t sort;
- Int_t tmp1, tmp2;
+ Int_t tmp1;
+ Int_t tmp2;
- if ( dir == 0 ) {
+ if (dir == 0) {
do {
sort = kTRUE;
- for ( Int_t i = 0; i < (nel-1); i++ )
- if ( x2[i+1] < x2[i] ) {
- tmp2 = x2[i];
- x2[i] = x2[i+1];
+ for (i = 0; i < (nel-1); i++) {
+ if (x2[i+1] < x2[i]) {
+ tmp2 = x2[i];
+ x2[i] = x2[i+1];
x2[i+1] = tmp2;
- tmp1 = x1[i];
- x1[i] = x1[i+1];
+ tmp1 = x1[i];
+ x1[i] = x1[i+1];
x1[i+1] = tmp1;
- sort = kFALSE;
+ sort = kFALSE;
}
- } while ( sort == kFALSE );
+ }
+ } while (sort == kFALSE);
}
- if ( dir == 1 ) {
+ if (dir == 1) {
do {
sort = kTRUE;
- for ( Int_t i = 0; i < (nel-1); i++ )
- if ( x2[i+1] > x2[i] ) {
- tmp2 = x2[i];
- x2[i] = x2[i+1];
+ for (i = 0; i < (nel-1); i++) {
+ if (x2[i+1] > x2[i]) {
+ tmp2 = x2[i];
+ x2[i] = x2[i+1];
x2[i+1] = tmp2;
- tmp1 = x1[i];
- x1[i] = x1[i+1];
+ tmp1 = x1[i];
+ x1[i] = x1[i+1];
x1[i+1] = tmp1;
- sort = kFALSE;
+ sort = kFALSE;
}
- } while ( sort == kFALSE );
+ }
+ } while (sort == kFALSE);
}
}
//_____________________________________________________________________________
-Bool_t AliTRDmcm::IsCluster(Float_t amp[3])
+Bool_t AliTRDmcm::IsCluster(Float_t amp[3]) const
{
//
// Find if the amplitudes amp[0], amp[1], amp[2] are a cluster
//
// -> shape
- if (amp[0] > amp[1] || amp[2] > amp[1]) return kFALSE;
+ if (amp[0] > amp[1] || amp[2] > amp[1]) {
+ return kFALSE;
+ }
// -> cluster amplitude
- if ((amp[0]+amp[1]+amp[2]) < fClusThr) return kFALSE;
+ if ((amp[0]+amp[1]+amp[2]) < fClusThr) {
+ return kFALSE;
+ }
// -> pad amplitude
- if (amp[0] < fPadThr && amp[2] < fPadThr) return kFALSE;
+ if (amp[0] < fPadThr && amp[2] < fPadThr) {
+ return kFALSE;
+ }
return kTRUE;
void AliTRDmcm::Filter(Int_t nexp, Int_t ftype)
{
//
- // exponential filter
+ // Exponential filter
//
+ Int_t iCol = 0;
+ Int_t iTime = 0;
+
Double_t sour[kMcmTBmax];
- Double_t Dtarg[kMcmTBmax];
- Int_t Itarg[kMcmTBmax];
+ Double_t dtarg[kMcmTBmax];
+ Int_t itarg[kMcmTBmax];
switch(ftype) {
- case 0:
+ case 0:
- for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
- for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
+ for (iCol = 0; iCol < kMcmCol; iCol++) {
+ for (iTime = 0; iTime < kMcmTBmax; iTime++) {
sour[iTime] = fADC[iCol][iTime];
}
- DeConvExpA(sour,Dtarg,kMcmTBmax,nexp);
- for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
- //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,Dtarg[iTime]);
- fADC[iCol][iTime] = TMath::Max(0.0,Dtarg[iTime]);
+ DeConvExpA(sour,dtarg,kMcmTBmax,nexp);
+ for (iTime = 0; iTime < kMcmTBmax; iTime++) {
+ fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
}
}
break;
case 1:
- for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
- for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
+ for (iCol = 0; iCol < kMcmCol; iCol++) {
+ for (iTime = 0; iTime < kMcmTBmax; iTime++) {
sour[iTime] = fADC[iCol][iTime];
}
- DeConvExpD(sour,Itarg,kMcmTBmax,nexp);
- for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
- fADC[iCol][iTime] = Itarg[iTime];
+ DeConvExpD(sour,itarg,kMcmTBmax,nexp);
+ for (iTime = 0; iTime < kMcmTBmax; iTime++) {
+ fADC[iCol][iTime] = itarg[iTime];
}
}
break;
case 2:
- for (Int_t iCol = 0; iCol < kMcmCol; iCol++) {
- for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
+ for (iCol = 0; iCol < kMcmCol; iCol++) {
+ for (iTime = 0; iTime < kMcmTBmax; iTime++) {
sour[iTime] = fADC[iCol][iTime];
}
- DeConvExpMI(sour,Dtarg,kMcmTBmax);
- for (Int_t iTime = 0; iTime < kMcmTBmax; iTime++) {
- //fADC[iCol][iTime] = (Int_t)TMath::Max(0.0,Dtarg[iTime]);
- fADC[iCol][iTime] = TMath::Max(0.0,Dtarg[iTime]);
+ DeConvExpMI(sour,dtarg,kMcmTBmax);
+ for (iTime = 0; iTime < kMcmTBmax; iTime++) {
+ fADC[iCol][iTime] = TMath::Max(0.0,dtarg[iTime]);
}
}
break;
default:
- printf("Invalid filter type %d ! \n",ftype);
+ AliError(Form("Invalid filter type %d ! \n",ftype));
return;
}
//_____________________________________________________________________________
void AliTRDmcm::DeConvExpA(Double_t *source, Double_t *target, Int_t n, Int_t nexp)
{
+ //
+ // Exponential filter "analog"
+ //
+ Int_t i = 0;
+ Int_t k = 0;
+ Double_t reminder[2];
+ Double_t correction;
+ Double_t result;
Double_t rates[2];
Double_t coefficients[2];
- // initialize (coefficient = alpha, rates = lambda)
+ // Initialize (coefficient = alpha, rates = lambda)
// FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
- Double_t R1, R2, C1, C2;
- R1 = (Double_t)fR1;
- R2 = (Double_t)fR2;
- C1 = (Double_t)fC1;
- C2 = (Double_t)fC2;
+ Double_t r1, r2, c1, c2;
+ r1 = (Double_t)fR1;
+ r2 = (Double_t)fR2;
+ c1 = (Double_t)fC1;
+ c2 = (Double_t)fC2;
- coefficients[0] = C1;
- coefficients[1] = C2;
+ coefficients[0] = c1;
+ coefficients[1] = c2;
- Double_t Dt = 0.100;
- rates[0] = TMath::Exp(-Dt/(R1));
- rates[1] = TMath::Exp(-Dt/(R2));
-
- Int_t i, k;
- Double_t reminder[2];
- Double_t correction, result;
+ Double_t dt = 0.1;
+ rates[0] = TMath::Exp(-dt/(r1));
+ rates[1] = TMath::Exp(-dt/(r2));
- /* attention: computation order is important */
- correction=0.0;
- for ( k=0; k<nexp; k++ ) reminder[k]=0.0;
+ // Attention: computation order is important
+ correction = 0.0;
+ for (k = 0; k < nexp; k++) {
+ reminder[k] = 0.0;
+ }
- for ( i=0; i<n; i++ ) {
- result = ( source[i] - correction ); // no rescaling
+ for (i = 0; i < n; i++) {
+
+ result = (source[i] - correction); // no rescaling
target[i] = result;
- for ( k=0; k<nexp; k++ ) reminder[k] = rates[k] * ( reminder[k] + coefficients[k] * result);
+ for (k = 0; k < nexp; k++) {
+ reminder[k] = rates[k] * (reminder[k] + coefficients[k] * result);
+ }
- correction=0.0;
- for ( k=0; k<nexp; k++ ) correction += reminder[k];
+ correction = 0.0;
+ for (k = 0; k < nexp; k++) {
+ correction += reminder[k];
+ }
+
}
}
//_____________________________________________________________________________
-void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp) {
-
- Int_t fAlpha_l, fAlpha_s, fLambda_l, fLambda_s, fTailPed;
- Int_t iAlpha_l, iAlpha_s, iLambda_l, iLambda_s;
+void AliTRDmcm::DeConvExpD(Double_t *source, Int_t *target, Int_t n, Int_t nexp)
+{
+ //
+ // Exponential filter "digital"
+ //
- // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
- // initialize (coefficient = alpha, rates = lambda)
+ Int_t i = 0;
- fLambda_l = 0; fAlpha_l = 0; fLambda_s = 0; fAlpha_s = 0;
- iLambda_l = 0; iAlpha_l = 0; iLambda_s = 0; iAlpha_s = 0;
+ Int_t fAlphaL;
+ Int_t fAlphaS;
+ Int_t fLambdaL;
+ Int_t fLambdaS;
+ Int_t fTailPed;
- Double_t Dt = 0.100;
+ Int_t iAlphaL;
+ Int_t iAlphaS;
+ Int_t iLambdaL;
+ Int_t iLambdaS;
- Double_t R1, R2, C1, C2;
- R1 = (Double_t)fR1;
- R2 = (Double_t)fR2;
- C1 = (Double_t)fC1;
- C2 = (Double_t)fC2;
+ // FilterOpt.C (aliroot@pel:/homel/aliroot/root/work/beamt/CERN02)
+ // initialize (coefficient = alpha, rates = lambda)
- fLambda_l = (Int_t)((TMath::Exp(-Dt/R1)-0.75)*2048.0);
- fLambda_s = (Int_t)((TMath::Exp(-Dt/R2)-0.25)*2048.0);
- iLambda_l = fLambda_l & 0x01FF; iLambda_l |= 0x0600; // 9 bit paramter + fixed bits
- iLambda_s = fLambda_s & 0x01FF; iLambda_s |= 0x0200; // 9 bit paramter + fixed bits
+ fLambdaL = 0;
+ fAlphaL = 0;
+ fLambdaS = 0;
+ fAlphaS = 0;
+ iLambdaL = 0;
+ iAlphaL = 0;
+ iLambdaS = 0;
+ iAlphaS = 0;
+
+ Double_t dt = 0.1;
+
+ Double_t r1;
+ Double_t r2;
+ Double_t c1;
+ Double_t c2;
+ r1 = (Double_t) fR1;
+ r2 = (Double_t) fR2;
+ c1 = (Double_t) fC1;
+ c2 = (Double_t) fC2;
+
+ fLambdaL = (Int_t)((TMath::Exp(-dt/r1) - 0.75) * 2048.0);
+ fLambdaS = (Int_t)((TMath::Exp(-dt/r2) - 0.25) * 2048.0);
+ iLambdaL = fLambdaL & 0x01FF; iLambdaL |= 0x0600; // 9 bit paramter + fixed bits
+ iLambdaS = fLambdaS & 0x01FF; iLambdaS |= 0x0200; // 9 bit paramter + fixed bits
if (nexp == 1) {
- fAlpha_l = (Int_t)(C1*2048.0);
- iAlpha_l = fAlpha_l & 0x03FF; // 10 bit paramter
+ fAlphaL = (Int_t) (c1 * 2048.0);
+ iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
}
if (nexp == 2) {
- fAlpha_l = (Int_t)(C1*2048.0);
- fAlpha_s = (Int_t)((C2-0.5)*2048.0);
- iAlpha_l = fAlpha_l & 0x03FF; // 10 bit paramter
- iAlpha_s = fAlpha_s & 0x03FF; iAlpha_s |= 0x0400; // 10 bit paramter + fixed bits
+ fAlphaL = (Int_t) (c1 * 2048.0);
+ fAlphaS = (Int_t) ((c2 - 0.5) * 2048.0);
+ iAlphaL = fAlphaL & 0x03FF; // 10 bit paramter
+ iAlphaS = fAlphaS & 0x03FF; iAlphaS |= 0x0400; // 10 bit paramter + fixed bits
}
- Double_t iAl = iAlpha_l / 2048.0; // alpha L: correspondence to floating point numbers
- Double_t iAs = iAlpha_s / 2048.0; // alpha S: correspondence to floating point numbers
- Double_t iLl = iLambda_l / 2048.0; // lambda L: correspondence to floating point numbers
- Double_t iLs = iLambda_s / 2048.0; // lambda S: correspondence to floating point numbers
-
- Int_t h1,h2;
- Int_t rem1, rem2;
- Int_t correction, result;
- Int_t iFactor = ((Int_t)fPedestal)<<2;
-
- Double_t xi = 1-(iLl*iAs + iLs*iAl); // calculation of equilibrium values of the
- rem1 = (Int_t)((iFactor/xi) * ((1-iLs)*iLl*iAl)); // internal registers to prevent switch on effects.
- rem2 = (Int_t)((iFactor/xi) * ((1-iLl)*iLs*iAs));
+ Double_t iAl = iAlphaL / 2048.0; // alpha L: correspondence to floating point numbers
+ Double_t iAs = iAlphaS / 2048.0; // alpha S: correspondence to floating point numbers
+ Double_t iLl = iLambdaL / 2048.0; // lambda L: correspondence to floating point numbers
+ Double_t iLs = iLambdaS / 2048.0; // lambda S: correspondence to floating point numbers
+
+ Int_t h1;
+ Int_t h2;
+ Int_t rem1;
+ Int_t rem2;
+ Int_t correction;
+ Int_t result;
+ Int_t iFactor = ((Int_t) fPedestal) << 2;
+
+ Double_t xi = 1 - (iLl*iAs + iLs*iAl); // Calculation of equilibrium values of the
+ rem1 = (Int_t) ((iFactor/xi) * ((1-iLs)*iLl*iAl)); // Internal registers to prevent switch on effects.
+ rem2 = (Int_t) ((iFactor/xi) * ((1-iLl)*iLs*iAs));
// further initialization
- if ( (rem1 + rem2) > 0x0FFF ) {
+ if ((rem1 + rem2) > 0x0FFF) {
correction = 0x0FFF;
- } else {
+ }
+ else {
correction = (rem1 + rem2) & 0x0FFF;
}
fTailPed = iFactor - correction;
- for (Int_t i=0; i<n; i++) {
+ for (i = 0; i < n; i++) {
- result = ( (Int_t)source[i] - correction );
- if ( result<0 ) {
+ result = ((Int_t)source[i] - correction);
+ if (result < 0) {
result = 0;
}
target[i] = result;
- h1 = ( rem1 + ((iAlpha_l * result) >> 11) );
- if ( h1 > 0x0FFF ) {
+ h1 = (rem1 + ((iAlphaL * result) >> 11));
+ if (h1 > 0x0FFF) {
h1 = 0x0FFF;
- } else {
+ }
+ else {
h1 &= 0x0FFF;
}
- h2 = ( rem2 + ((iAlpha_s * result) >> 11));
- if ( h2 > 0x0FFF ) {
+ h2 = (rem2 + ((iAlphaS * result) >> 11));
+ if (h2 > 0x0FFF) {
h2 = 0x0FFF;
- } else {
+ }
+ else {
h2 &= 0x0FFF;
}
- rem1 = (iLambda_l * h1 ) >> 11;
- rem2 = (iLambda_s * h2 ) >> 11;
+ rem1 = (iLambdaL * h1 ) >> 11;
+ rem2 = (iLambdaS * h2 ) >> 11;
- if ( (rem1 + rem2) > 0x0FFF ) {
+ if ((rem1 + rem2) > 0x0FFF) {
correction = 0x0FFF;
- } else {
+ }
+ else {
correction = (rem1 + rem2) & 0x0FFF;
}
+
}
}
//_____________________________________________________________________________
-void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n) {
+void AliTRDmcm::DeConvExpMI(Double_t *source, Double_t *target, Int_t n)
+{
+ //
+ // Exponential filter (M. Ivanov)
+ //
- Double_t Sig1[100], Sig2[100], Sig3[100];//, Sig4[100];
- for (Int_t i = 0; i < n; i++) Sig1[i] = source[i];
+ Int_t i = 0;
- Float_t Dt = 0.100;
+ Double_t sig1[100];
+ Double_t sig2[100];
+ Double_t sig3[100];
- //Float_t lambda0 = 9.8016*Dt; // short
- //Float_t lambda1 = 1.0778*Dt; // long
+ for (i = 0; i < n; i++) {
+ sig1[i] = source[i];
+ }
- Float_t lambda0 = (1.0/fR2)*Dt;
- Float_t lambda1 = (1.0/fR1)*Dt;
+ Float_t dt = 0.1;
- TailMakerSpline(Sig1,Sig2,lambda0,n);
- TailCancelationMI(Sig2,Sig3,0.7,lambda1,n);
+ Float_t lambda0 = (1.0 / fR2) * dt;
+ Float_t lambda1 = (1.0 / fR1) * dt;
- for (Int_t i = 0; i < n; i++) target[i] = Sig3[i];
+ TailMakerSpline(sig1,sig2,lambda0,n);
+ TailCancelationMI(sig2,sig3,0.7,lambda1,n);
+
+ for (i = 0; i < n; i++) {
+ target[i] = sig3[i];
+ }
}
//______________________________________________________________________________
-void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n) {
-
- Double_t l = TMath::Exp(-lambda*0.5);
- //
- //
- Double_t in[1000];
- Double_t out[1000];
- // initialize in[] and out[] goes 0 ... 2*n+19
- for (Int_t i=0; i<n*2+20; i++){
- in[i] = 0;
- out[i]= 0;
- }
-
- // in[] goes 0, 1
- in[0] = ampin[0];
- in[1] = (ampin[0]+ampin[1])*0.5;
+void AliTRDmcm::TailMakerSpline(Double_t *ampin, Double_t *ampout, Double_t lambda, Int_t n)
+{
+ //
+ // Special filter (M. Ivanov)
+ //
+
+ Int_t i = 0;
+
+ Double_t l = TMath::Exp(-lambda*0.5);
+ Double_t in[1000];
+ Double_t out[1000];
+
+ // Initialize in[] and out[] goes 0 ... 2*n+19
+ for (i = 0; i < n*2+20; i++) {
+ in[i] = 0;
+ out[i] = 0;
+ }
+
+ // in[] goes 0, 1
+ in[0] = ampin[0];
+ in[1] = (ampin[0] + ampin[1]) * 0.5;
- // add charge to the end
- for (Int_t i=0; i<22; i++){
- // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
- in[2*(n-1)+i] = ampin[n-1];
- }
-
- // use arithm mean
- for (Int_t i=1; i<n-1; i++){
- // in[] goes 2, 3, ... , 2*n-4, 2*n-3
- in[2*i] = ampin[i];
- in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
- }
- /*
- // add charge to the end
- for (Int_t i=-2; i<22; i++){
- // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
- in[2*(n-1)+i] = ampin[n-1];
- }
-
- // use spline mean
- for (Int_t i=1; i<n-2; i++){
- // in[] goes 2, 3, ... , 2*n-6, 2*n-5
- in[2*i] = ampin[i];
- in[2*i+1] = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16;
- }
- */
- //
- Double_t temp;
- out[2*n] = in[2*n];
- temp = 0;
- for (int i=2*n; i>=0; i--){
- out[i] = in[i] + temp;
- temp = l*(temp+in[i]);
- }
-
- //
- for (int i=0;i<n;i++){
- //ampout[i] = out[2*i+1]; // org
- ampout[i] = out[2*i];
- }
+ // Add charge to the end
+ for (i = 0; i < 22; i++) {
+ // in[] goes 2*n-2, 2*n-1, ... , 2*n+19
+ in[2*(n-1)+i] = ampin[n-1];
+ }
+
+ // Use arithmetic mean
+ for (i = 1; i < n-1; i++) {
+ // in[] goes 2, 3, ... , 2*n-4, 2*n-3
+ in[2*i] = ampin[i];
+ in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
+ }
+
+ Double_t temp;
+ out[2*n] = in[2*n];
+ temp = 0;
+ for (i = 2*n; i >= 0; i--) {
+ out[i] = in[i] + temp;
+ temp = l*(temp+in[i]);
+ }
+
+ for (i = 0; i < n; i++){
+ //ampout[i] = out[2*i+1]; // org
+ ampout[i] = out[2*i];
+ }
}
//______________________________________________________________________________
-void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout, Double_t norm, Double_t lambda, Int_t n) {
-
- Double_t L = TMath::Exp(-lambda*0.5);
- Double_t K = L*(1.-norm*lambda*0.5);
- //
- //
- Double_t in[1000];
- Double_t out[1000];
- // initialize in[] and out[] goes 0 ... 2*n+19
- for (Int_t i=0; i<n*2+20; i++){
- in[i] = 0;
- out[i]= 0;
- }
-
- // in[] goes 0, 1
- in[0] = ampin[0];
- in[1] = (ampin[0]+ampin[1])*0.5;
-
- // add charge to the end
- for (Int_t i=-2; i<22; i++){
- // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
- in[2*(n-1)+i] = ampin[n-1];
- }
-
- //
- for (Int_t i=1; i<n-2; i++){
- // in[] goes 2, 3, ... , 2*n-6, 2*n-5
- in[2*i] = ampin[i];
- in[2*i+1] = (9.*(ampin[i]+ampin[i+1])-(ampin[i-1]+ampin[i+2]))/16.;
- //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.;
- }
-
- //
- Double_t temp;
- out[0] = in[0];
- temp = in[0];
- for (int i=1; i<=2*n; i++){
- out[i] = in[i] + (K-L)*temp;
- temp = in[i] + K*temp;
- }
- //
- //
- for (int i=0; i<n; i++){
- //ampout[i] = out[2*i+1]; // org
- //ampout[i] = TMath::Max(out[2*i+1], 0.0); // org
- ampout[i] = TMath::Max(out[2*i], 0.0);
- }
+void AliTRDmcm::TailCancelationMI(Double_t *ampin, Double_t *ampout
+ , Double_t norm, Double_t lambda, Int_t n)
+{
+ //
+ // Special filter (M. Ivanov)
+ //
+
+ Int_t i = 0;
+
+ Double_t l = TMath::Exp(-lambda*0.5);
+ Double_t k = l*(1.0 - norm*lambda*0.5);
+ Double_t in[1000];
+ Double_t out[1000];
+
+ // Initialize in[] and out[] goes 0 ... 2*n+19
+ for (i = 0; i < n*2+20; i++) {
+ in[i] = 0;
+ out[i] = 0;
+ }
+
+ // in[] goes 0, 1
+ in[0] = ampin[0];
+ in[1] = (ampin[0]+ampin[1])*0.5;
+
+ // Add charge to the end
+ for (i =-2; i < 22; i++) {
+ // in[] goes 2*n-4, 2*n-3, ... , 2*n+19
+ in[2*(n-1)+i] = ampin[n-1];
+ }
+
+ for (i = 1; i < n-2; i++) {
+ // in[] goes 2, 3, ... , 2*n-6, 2*n-5
+ in[2*i] = ampin[i];
+ in[2*i+1] = (9.0 * (ampin[i]+ampin[i+1]) - (ampin[i-1]+ampin[i+2])) / 16.0;
+ //in[2*i+1] = ((ampin[i]+ampin[i+1]))/2.0;
+ }
+
+ Double_t temp;
+ out[0] = in[0];
+ temp = in[0];
+ for (i = 1; i <= 2*n; i++) {
+ out[i] = in[i] + (k-l)*temp;
+ temp = in[i] + k *temp;
+ }
+
+ for (i = 0; i < n; i++) {
+ //ampout[i] = out[2*i+1]; // org
+ //ampout[i] = TMath::Max(out[2*i+1],0.0); // org
+ ampout[i] = TMath::Max(out[2*i],0.0);
+ }
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff