X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=TRD%2FAliTRDclusterizerV1.cxx;h=3e19c295ccbedba0e0c482886c0d52527cb0a048;hb=9679d5dbb99896e8003cf8573df706e40e172176;hp=72ba361286aceea7d0b440a24e81da653838699b;hpb=6f1e466d6db177aa095cc30f57958d5e62389e1c;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDclusterizerV1.cxx b/TRD/AliTRDclusterizerV1.cxx index 72ba361286a..3e19c295ccb 100644 --- a/TRD/AliTRDclusterizerV1.cxx +++ b/TRD/AliTRDclusterizerV1.cxx @@ -1,3 +1,4 @@ + /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -13,15 +14,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.1.4.1 2000/05/08 15:09:01 cblume -Introduce AliTRDdigitsManager - -Revision 1.1 2000/02/28 18:58:54 cblume -Add new TRD classes - -*/ +/* $Id$ */ /////////////////////////////////////////////////////////////////////////////// // // @@ -30,13 +23,28 @@ Add new TRD classes /////////////////////////////////////////////////////////////////////////////// #include +#include +#include +#include + +#include "AliRun.h" +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliRawReader.h" #include "AliTRDclusterizerV1.h" #include "AliTRDmatrix.h" #include "AliTRDgeometry.h" -#include "AliTRDdigitizer.h" -#include "AliTRDrecPoint.h" #include "AliTRDdataArrayF.h" +#include "AliTRDdataArrayI.h" +#include "AliTRDdigitsManager.h" +#include "AliTRDpadPlane.h" +#include "AliTRDrawData.h" +#include "AliTRDcalibDB.h" +#include "AliTRDSimParam.h" +#include "AliTRDRecParam.h" +#include "AliTRDCommonParam.h" +#include "AliTRDcluster.h" ClassImp(AliTRDclusterizerV1) @@ -47,7 +55,7 @@ AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer() // AliTRDclusterizerV1 default constructor // - fDigitsManager = NULL; + fDigitsManager = 0; } @@ -60,32 +68,58 @@ AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title // fDigitsManager = new AliTRDdigitsManager(); + fDigitsManager->CreateArrays(); + +} + +//_____________________________________________________________________________ +AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c) +:AliTRDclusterizer(c) +{ + // + // AliTRDclusterizerV1 copy constructor + // - Init(); + ((AliTRDclusterizerV1 &) c).Copy(*this); } //_____________________________________________________________________________ AliTRDclusterizerV1::~AliTRDclusterizerV1() { + // + // AliTRDclusterizerV1 destructor + // if (fDigitsManager) { delete fDigitsManager; + fDigitsManager = NULL; } } //_____________________________________________________________________________ -void AliTRDclusterizerV1::Init() +AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c) +{ + // + // Assignment operator + // + + if (this != &c) ((AliTRDclusterizerV1 &) c).Copy(*this); + return *this; + +} + +//_____________________________________________________________________________ +void AliTRDclusterizerV1::Copy(TObject &c) const { // - // Initializes the cluster finder + // Copy function // - // The default parameter for the clustering - fClusMaxThresh = 5.0; - fClusSigThresh = 2.0; - fClusMethod = 1; + ((AliTRDclusterizerV1 &) c).fDigitsManager = 0; + + AliTRDclusterizer::Copy(c); } @@ -96,19 +130,37 @@ Bool_t AliTRDclusterizerV1::ReadDigits() // Reads the digits arrays from the input aliroot file // - if (!fInputFile) { - printf("AliTRDclusterizerV1::ReadDigits -- "); + if (!fRunLoader) { + printf(" "); printf("No input file open\n"); return kFALSE; } + AliLoader* loader = fRunLoader->GetLoader("TRDLoader"); + if (!loader->TreeD()) loader->LoadDigits(); // Read in the digit arrays - return (fDigitsManager->ReadDigits()); + return (fDigitsManager->ReadDigits(loader->TreeD())); } //_____________________________________________________________________________ -Bool_t AliTRDclusterizerV1::MakeCluster() +Bool_t AliTRDclusterizerV1::ReadDigits(AliRawReader* rawReader) +{ + // + // Reads the digits arrays from the ddl file + // + + AliTRDrawData raw; + raw.SetDebug(1); + + fDigitsManager = raw.Raw2Digits(rawReader); + + return kTRUE; + +} + +//_____________________________________________________________________________ +Bool_t AliTRDclusterizerV1::MakeClusters() { // // Generates the cluster. @@ -116,236 +168,389 @@ Bool_t AliTRDclusterizerV1::MakeCluster() Int_t row, col, time; - // Get the pointer to the detector class and check for version 1 - AliTRD *TRD = (AliTRD*) gAlice->GetDetector("TRD"); - if (TRD->IsVersion() != 1) { - printf("AliTRDclusterizerV1::MakeCluster -- "); + /* + if (fTRD->IsVersion() != 1) { + printf(" "); printf("TRD must be version 1 (slow simulator).\n"); return kFALSE; } + */ // Get the geometry - AliTRDgeometry *Geo = TRD->GetGeometry(); - - printf("AliTRDclusterizerV1::MakeCluster -- "); - printf("Start creating clusters.\n"); - - AliTRDdataArrayI *Digits; - - // Parameters - Float_t maxThresh = fClusMaxThresh; // threshold value for maximum - Float_t signalThresh = fClusSigThresh; // threshold value for digit signal - Int_t clusteringMethod = fClusMethod; // clustering method option (for testing) - + AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader); + AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); + if (!calibration) + { + printf(" "); + printf("ERROR getting instance of AliTRDcalibDB"); + return kFALSE; + } + + AliTRDSimParam* simParam = AliTRDSimParam::Instance(); + if (!simParam) + { + printf(" "); + printf("ERROR getting instance of AliTRDSimParam"); + return kFALSE; + } + + AliTRDRecParam* recParam = AliTRDRecParam::Instance(); + if (!recParam) + { + printf(" "); + printf("ERROR getting instance of AliTRDRecParam"); + return kFALSE; + } + + AliTRDCommonParam* commonParam = AliTRDCommonParam::Instance(); + if (!commonParam) + { + printf(" "); + printf("Could not get common params\n"); + return kFALSE; + } + + Float_t ADCthreshold = simParam->GetADCthreshold(); + + if (fVerbose > 0) { + //printf(" "); + //printf("OmegaTau = %f \n",omegaTau); + printf(" "); + printf("Start creating clusters.\n"); + } + + AliTRDdataArrayI *digitsIn; + AliTRDdataArrayI *track0; + AliTRDdataArrayI *track1; + AliTRDdataArrayI *track2; + + // Threshold value for the maximum + Float_t maxThresh = recParam->GetClusMaxThresh(); + // Threshold value for the digit signal + Float_t sigThresh = recParam->GetClusSigThresh(); // Iteration limit for unfolding procedure - const Float_t epsilon = 0.01; + const Float_t kEpsilon = 0.01; - const Int_t nClus = 3; - const Int_t nSig = 5; + const Int_t kNclus = 3; + const Int_t kNsig = 5; + const Int_t kNtrack = 3 * kNclus; - Int_t chamBeg = 0; - Int_t chamEnd = kNcham; - if (TRD->GetSensChamber() >= 0) { - chamBeg = TRD->GetSensChamber(); - chamEnd = chamBeg + 1; - } - Int_t planBeg = 0; - Int_t planEnd = kNplan; - if (TRD->GetSensPlane() >= 0) { - planBeg = TRD->GetSensPlane(); - planEnd = planBeg + 1; - } - Int_t sectBeg = 0; - Int_t sectEnd = kNsect; - if (TRD->GetSensSector() >= 0) { - sectBeg = TRD->GetSensSector(); - sectEnd = sectBeg + 1; + Int_t iType = 0; + Int_t iUnfold = 0; + Double_t ratioLeft = 1.0; + Double_t ratioRight = 1.0; + + // + Double_t padSignal[kNsig]; + Double_t clusterSignal[kNclus]; + Double_t clusterPads[kNclus]; + Int_t clusterTracks[kNtrack]; + + Int_t chamBeg = 0; + Int_t chamEnd = AliTRDgeometry::Ncham(); + Int_t planBeg = 0; + Int_t planEnd = AliTRDgeometry::Nplan(); + Int_t sectBeg = 0; + Int_t sectEnd = AliTRDgeometry::Nsect(); + + Int_t nTimeTotal = calibration->GetNumberOfTimeBins(); + + if (fVerbose > 0) { + printf(" "); + printf("Number of Time Bins = %d.\n",nTimeTotal); } - // *** Start clustering *** in every chamber + // Start clustering in every chamber for (Int_t icham = chamBeg; icham < chamEnd; icham++) { for (Int_t iplan = planBeg; iplan < planEnd; iplan++) { for (Int_t isect = sectBeg; isect < sectEnd; isect++) { - Int_t idet = Geo->GetDetector(iplan,icham,isect); - - Int_t nClusters = 0; - printf("AliTRDclusterizerV1::MakeCluster -- "); - printf("Analyzing chamber %d, plane %d, sector %d.\n" - ,icham,iplan,isect); - - Int_t nRowMax = Geo->GetRowMax(iplan,icham,isect); - Int_t nColMax = Geo->GetColMax(iplan); - Int_t nTimeMax = Geo->GetTimeMax(); - - // Create a detector matrix to keep maxima - AliTRDmatrix *digitMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax - ,isect,icham,iplan); - // Create a matrix to contain maximum flags - AliTRDmatrix *maximaMatrix = new AliTRDmatrix(nRowMax,nColMax,nTimeMax - ,isect,icham,iplan); - - // Read in the digits - Digits = fDigitsManager->GetDigits(idet); - - // Loop through the detector pixel - for (time = 0; time < nTimeMax; time++) { - for ( col = 0; col < nColMax; col++) { - for ( row = 0; row < nRowMax; row++) { - - Int_t signal = Digits->GetData(row,col,time); - Int_t index = Digits->GetIndex(row,col,time); - - // Fill the detector matrix - if (signal > signalThresh) { - // Store the signal amplitude - digitMatrix->SetSignal(row,col,time,signal); - // Store the digits number - digitMatrix->AddTrack(row,col,time,index); - } + Int_t idet = geo->GetDetector(iplan,icham,isect); - } - } - } - - // Loop chamber and find maxima in digitMatrix - for ( row = 0; row < nRowMax; row++) { - for ( col = 1; col < nColMax; col++) { - for (time = 0; time < nTimeMax; time++) { - - if (digitMatrix->GetSignal(row,col,time) - < digitMatrix->GetSignal(row,col - 1,time)) { - // really maximum? - if (col > 1) { - if (digitMatrix->GetSignal(row,col - 2,time) - < digitMatrix->GetSignal(row,col - 1,time)) { - // yes, so set maximum flag - maximaMatrix->SetSignal(row,col - 1,time,1); - } - else maximaMatrix->SetSignal(row,col - 1,time,0); - } - } - - } // time - } // col - } // row - - // now check maxima and calculate cluster position - for ( row = 0; row < nRowMax; row++) { - for ( col = 1; col < nColMax; col++) { - for (time = 0; time < nTimeMax; time++) { - - if ((maximaMatrix->GetSignal(row,col,time) > 0) - && (digitMatrix->GetSignal(row,col,time) > maxThresh)) { - - // Ratio resulting from unfolding - Float_t ratio = 0; - // Signals on max and neighbouring pads - Float_t padSignal[nSig] = {0}; - // Signals from cluster - Float_t clusterSignal[nClus] = {0}; - // Cluster pad info - Float_t clusterPads[nClus] = {0}; - // Cluster digit info - Int_t clusterDigit[nClus] = {0}; - - for (Int_t iPad = 0; iPad < nClus; iPad++) { - clusterSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time); - clusterDigit[iPad] = digitMatrix->GetTrack(row,col-1+iPad,time,0); - } - - // neighbouring maximum on right side? - if (col < nColMax - 2) { - if (maximaMatrix->GetSignal(row,col + 2,time) > 0) { + Int_t nRowMax = commonParam->GetRowMax(iplan,icham,isect); + Int_t nColMax = commonParam->GetColMax(iplan); - for (Int_t iPad = 0; iPad < 5; iPad++) { - padSignal[iPad] = digitMatrix->GetSignal(row,col-1+iPad,time); - } + Int_t nClusters = 0; + Int_t nClusters2pad = 0; + Int_t nClusters3pad = 0; + Int_t nClusters4pad = 0; + Int_t nClusters5pad = 0; + Int_t nClustersLarge = 0; - // unfold: - ratio = Unfold(epsilon, padSignal); - - // set signal on overlapping pad to ratio - clusterSignal[2] *= ratio; + if (fVerbose > 0) { + printf(" "); + printf("Analyzing chamber %d, plane %d, sector %d.\n" + ,icham,iplan,isect); + } - } + AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham); + + // Get the digits + digitsIn = fDigitsManager->GetDigits(idet); + digitsIn->Expand(); + AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow(), digitsIn->GetNcol(), digitsIn->GetNtime()); + + Transform(digitsIn, digitsOut, idet, nRowMax, nColMax, nTimeTotal, ADCthreshold); + + track0 = fDigitsManager->GetDictionary(idet,0); + track0->Expand(); + track1 = fDigitsManager->GetDictionary(idet,1); + track1->Expand(); + track2 = fDigitsManager->GetDictionary(idet,2); + track2->Expand(); + + // Loop through the chamber and find the maxima + for ( row = 0; row < nRowMax; row++) { + for ( col = 2; col < nColMax; col++) { + //for ( col = 4; col < nColMax-2; col++) { + for (time = 0; time < nTimeTotal; time++) { + + Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); + Float_t signalM = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); + Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)); + +// // Lonok for the maximum +// if (signalM >= maxThresh) { +// if (((signalL >= sigThresh) && +// (signalL < signalM)) || +// ((signalR >= sigThresh) && +// (signalR < signalM))) { +// // Maximum found, mark the position by a negative signal +// digitsOut->SetDataUnchecked(row,col-1,time,-signalM); +// } +// } + // Look for the maximum + if (signalM >= maxThresh) { + if ( (TMath::Abs(signalL)<=signalM) && (TMath::Abs(signalR)<=signalM) && + (TMath::Abs(signalL)+TMath::Abs(signalR))>sigThresh ) { + // Maximum found, mark the position by a negative signal + digitsOut->SetDataUnchecked(row,col-1,time,-signalM); + } + } + } + } + } + + // Now check the maxima and calculate the cluster position + for ( row = 0; row < nRowMax ; row++) { + for (time = 0; time < nTimeTotal; time++) { + for ( col = 1; col < nColMax-1; col++) { + + // Maximum found ? + if (digitsOut->GetDataUnchecked(row,col,time) < 0) { + + Int_t iPad; + for (iPad = 0; iPad < kNclus; iPad++) { + Int_t iPadCol = col - 1 + iPad; + clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row + ,iPadCol + ,time)); + clusterTracks[3*iPad ] = track0->GetDataUnchecked(row,iPadCol,time) - 1; + clusterTracks[3*iPad+1] = track1->GetDataUnchecked(row,iPadCol,time) - 1; + clusterTracks[3*iPad+2] = track2->GetDataUnchecked(row,iPadCol,time) - 1; } - - // Calculate the position of the cluster - switch (clusteringMethod) { - case 1: - // method 1: simply center of mass - clusterPads[0] = row + 0.5; - clusterPads[1] = col - 0.5 + (clusterSignal[2] - clusterSignal[0]) / - (clusterSignal[0] + clusterSignal[1] + clusterSignal[2]); - clusterPads[2] = time + 0.5; - - nClusters++; - break; - case 2: - // method 2: integral gauss fit on 3 pads - TH1F *hPadCharges = new TH1F("hPadCharges", "Charges on center 3 pads" - , 5, -1.5, 3.5); - for (Int_t iCol = -1; iCol <= 3; iCol++) { - if (clusterSignal[iCol] < 1) clusterSignal[iCol] = 1; - hPadCharges->Fill(iCol, clusterSignal[iCol]); - } - hPadCharges->Fit("gaus", "IQ", "SAME", -0.5, 2.5); - TF1 *fPadChargeFit = hPadCharges->GetFunction("gaus"); - Double_t colMean = fPadChargeFit->GetParameter(1); - clusterPads[0] = row + 0.5; - clusterPads[1] = col - 1.5 + colMean; - clusterPads[2] = time + 0.5; - - delete hPadCharges; - - nClusters++; + // Count the number of pads in the cluster + Int_t nPadCount = 0; + Int_t ii = 0; + while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) + >= sigThresh) { + nPadCount++; + ii++; + if (col-ii < 0) break; + } + ii = 0; + while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) + >= sigThresh) { + nPadCount++; + ii++; + if (col+ii+1 >= nColMax) break; + } + + nClusters++; + switch (nPadCount) { + case 2: + iType = 0; + nClusters2pad++; + break; + case 3: + iType = 1; + nClusters3pad++; break; + case 4: + iType = 2; + nClusters4pad++; + break; + case 5: + iType = 3; + nClusters5pad++; + break; + default: + iType = 4; + nClustersLarge++; + break; + }; + + // Look for 5 pad cluster with minimum in the middle + Bool_t fivePadCluster = kFALSE; + if (col < nColMax-3) { + if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) { + fivePadCluster = kTRUE; + } + if ((fivePadCluster) && (col < nColMax-5)) { + if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) { + fivePadCluster = kFALSE; + } + } + if ((fivePadCluster) && (col > 1)) { + if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) { + fivePadCluster = kFALSE; + } + } + } + + // 5 pad cluster + // Modify the signal of the overlapping pad for the left part + // of the cluster which remains from a previous unfolding + if (iUnfold) { + clusterSignal[0] *= ratioLeft; + iType = 5; + iUnfold = 0; + } + + // Unfold the 5 pad cluster + if (fivePadCluster) { + for (iPad = 0; iPad < kNsig; iPad++) { + padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row + ,col-1+iPad + ,time)); + } + // Unfold the two maxima and set the signal on + // the overlapping pad to the ratio + ratioRight = Unfold(kEpsilon,iplan,padSignal); + ratioLeft = 1.0 - ratioRight; + clusterSignal[2] *= ratioRight; + iType = 5; + iUnfold = 1; } - Float_t clusterCharge = clusterSignal[0] - + clusterSignal[1] - + clusterSignal[2]; + Double_t clusterCharge = clusterSignal[0] + + clusterSignal[1] + + clusterSignal[2]; + + // The position of the cluster + clusterPads[0] = row + 0.5; + // Take the shift of the additional time bins into account + clusterPads[2] = time + 0.5; + + if (recParam->LUTOn()) { + // Calculate the position of the cluster by using the + // lookup table method + clusterPads[1] = recParam->LUTposition(iplan,clusterSignal[0] + ,clusterSignal[1] + ,clusterSignal[2]); + } + else { + // Calculate the position of the cluster by using the + // center of gravity method + for (Int_t i=0;i<5;i++) padSignal[i]=0; + padSignal[2] = TMath::Abs(digitsOut->GetDataUnchecked(row,col,time)); // central pad + padSignal[1] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-1,time)); // left pad + padSignal[3] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+1,time)); // right pad + if (col>2 &&TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)GetDataUnchecked(row,col-2,time)); + } + if (colGetDataUnchecked(row,col+2,time)GetDataUnchecked(row,col+2,time)); + } + clusterPads[1] = GetCOG(padSignal); + + } + + Double_t q0 = clusterSignal[0]; + Double_t q1 = clusterSignal[1]; + Double_t q2 = clusterSignal[2]; + Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) / + (clusterCharge*clusterCharge); + + + + // Calculate the position and the error + + // correct for t0 + Int_t clusterTimeBin = TMath::Nint(time - calibration->GetT0(idet, col, row)); + + Double_t colSize = padPlane->GetColSize(col); + Double_t rowSize = padPlane->GetRowSize(row); + Double_t clusterPos[3]; + clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1]+0.5)*colSize; // MI change + clusterPos[1] = padPlane->GetRowPos(row) - 0.5*rowSize; //MI change + clusterPos[2] = CalcXposFromTimebin(clusterPads[2], idet, col, row); + Double_t clusterSig[2]; + clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize; + clusterSig[1] = rowSize * rowSize / 12.; + + // Add the cluster to the output array - TRD->AddRecPoint(clusterPads,clusterDigit,idet,clusterCharge); - + AliTRDcluster * cluster = AddCluster(clusterPos + ,clusterTimeBin + ,idet + ,clusterCharge + ,clusterTracks + ,clusterSig + ,iType,clusterPads[1]); + // + // + Short_t signals[7]={0,0,0,0,0,0,0}; + for (Int_t jPad = col-3;jPad<=col+3;jPad++){ + if (jPad<0 ||jPad>=nColMax-1) continue; + signals[jPad-col+3] = TMath::Nint(TMath::Abs(digitsOut->GetDataUnchecked(row,jPad,time))); + } + cluster->SetSignals(signals); } - } // time - } // col - } // row - - printf("AliTRDclusterizerV1::MakeCluster -- "); - printf("Number of clusters found: %d\n",nClusters); - - delete digitMatrix; - delete maximaMatrix; - - } // isect - } // iplan - } // icham - - printf("AliTRDclusterizerV1::MakeCluster -- "); - printf("Total number of points found: %d\n" - ,TRD->RecPoints()->GetEntries()); - - // Get the pointer to the cluster branch - TTree *ClusterTree = gAlice->TreeR(); + } + } + } + + delete digitsOut; + + // Compress the arrays + track0->Compress(1,0); + track1->Compress(1,0); + track2->Compress(1,0); + + // Write the cluster and reset the array + WriteClusters(idet); + ResetRecPoints(); + } + } + } + + if (fVerbose > 0) { + printf(" "); + printf("Done.\n"); + } - // Fill the cluster-branch - printf("AliTRDclusterizerV1::MakeCluster -- "); - printf("Fill the cluster tree.\n"); - ClusterTree->Fill(); - printf("AliTRDclusterizerV1::MakeCluster -- "); - printf("Done.\n"); + //delete digitsIn; return kTRUE; } //_____________________________________________________________________________ -Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal) +Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5]) +{ + // + // get COG position + // used for clusters with more than 3 pads - where LUT not applicable + Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4]; + Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum; + return res; +} + +//_____________________________________________________________________________ +Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal) { // // Method to unfold neighbouring maxima. @@ -354,41 +559,49 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal) // The resulting ratio is then returned to the calling method. // - Int_t itStep = 0; // count iteration steps + AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); + if (!calibration) + { + printf(" "); + printf("ERROR getting instance of AliTRDcalibDB"); + return kFALSE; + } + + Int_t irc = 0; + Int_t itStep = 0; // Count iteration steps - Float_t ratio = 0.5; // start value for ratio - Float_t prevRatio = 0; // store previous ratio + Double_t ratio = 0.5; // Start value for ratio + Double_t prevRatio = 0; // Store previous ratio - Float_t newLeftSignal[3] = {0}; // array to store left cluster signal - Float_t newRightSignal[3] = {0}; // array to store right cluster signal + Double_t newLeftSignal[3] = {0}; // Array to store left cluster signal + Double_t newRightSignal[3] = {0}; // Array to store right cluster signal + Double_t newSignal[3] = {0}; - // start iteration: + // Start the iteration while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { itStep++; prevRatio = ratio; - // cluster position according to charge ratio - Float_t maxLeft = (ratio*padSignal[2] - padSignal[0]) / - (padSignal[0] + padSignal[1] + ratio*padSignal[2]); - Float_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) / - ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]); - - // set cluster charge ratio - Float_t ampLeft = padSignal[1]; - Float_t ampRight = padSignal[3]; + // Cluster position according to charge ratio + Double_t maxLeft = (ratio*padSignal[2] - padSignal[0]) + / (padSignal[0] + padSignal[1] + ratio*padSignal[2]); + Double_t maxRight = (padSignal[4] - (1-ratio)*padSignal[2]) + / ((1-ratio)*padSignal[2] + padSignal[3] + padSignal[4]); - // apply pad response to parameters - newLeftSignal[0] = ampLeft*PadResponse(-1 - maxLeft); - newLeftSignal[1] = ampLeft*PadResponse( 0 - maxLeft); - newLeftSignal[2] = ampLeft*PadResponse( 1 - maxLeft); + // Set cluster charge ratio + irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal); + Double_t ampLeft = padSignal[1] / newSignal[1]; + irc = calibration->PadResponse(1.0,maxRight,plane,newSignal); + Double_t ampRight = padSignal[3] / newSignal[1]; - newRightSignal[0] = ampRight*PadResponse(-1 - maxRight); - newRightSignal[1] = ampRight*PadResponse( 0 - maxRight); - newRightSignal[2] = ampRight*PadResponse( 1 - maxRight); + // Apply pad response to parameters + irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal ); + irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal); - // calculate new overlapping ratio - ratio = newLeftSignal[2]/(newLeftSignal[2] + newRightSignal[0]); + // Calculate new overlapping ratio + ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] / + (newLeftSignal[2] + newRightSignal[0])); } @@ -397,22 +610,140 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Float_t* padSignal) } //_____________________________________________________________________________ -Float_t AliTRDclusterizerV1::PadResponse(Float_t x) +void AliTRDclusterizerV1::Transform(AliTRDdataArrayI* digitsIn, + AliTRDdataArrayF* digitsOut, + Int_t idet, Int_t nRowMax, + Int_t nColMax, Int_t nTimeTotal, + Float_t ADCthreshold) +{ + + // + // Apply gain factor + // Apply tail cancellation: Transform digitsIn to digitsOut + // + + + AliTRDRecParam* recParam = AliTRDRecParam::Instance(); + if (!recParam) + { + printf(" "); + printf("ERROR getting instance of AliTRDRecParam"); + return; + } + AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); + + Double_t *inADC = new Double_t[nTimeTotal]; // adc data before tail cancellation + Double_t *outADC = new Double_t[nTimeTotal]; // adc data after tail cancellation + + if (fVerbose > 0) { + printf(" "); + printf("Tail cancellation (nExp = %d) for detector %d.\n", + recParam->GetTCnexp(),idet); + } + + for (Int_t iRow = 0; iRow < nRowMax; iRow++ ) { + for (Int_t iCol = 0; iCol < nColMax; iCol++ ) { + for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) { + // + // add gain + // + Double_t gain = calibration->GetGainFactor(idet, iCol, iRow); + if (gain==0) { + AliError("Not a valid gain\n"); + } + inADC[iTime] = digitsIn->GetDataUnchecked(iRow, iCol, iTime); + + inADC[iTime] /= gain; + outADC[iTime] = inADC[iTime]; + } + + // Apply the tail cancelation via the digital filter + if (recParam->TCOn()) + { + DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp()); + } + + for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) { + // Store the amplitude of the digit if above threshold + if (outADC[iTime] > ADCthreshold) { + if (fVerbose > 1) + { + printf(" iRow = %d, iCol = %d, iTime = %d, adc = %f\n" + ,iRow,iCol,iTime,outADC[iTime]); + } + digitsOut->SetDataUnchecked(iRow,iCol,iTime,outADC[iTime]); + } + + } + + } + + } + + delete [] inADC; + delete [] outADC; + + return; + +} + + +//_____________________________________________________________________________ +void AliTRDclusterizerV1::DeConvExp(Double_t *source, Double_t *target, + Int_t n, Int_t nexp) { // - // The pad response for the chevron pads. - // We use a simple Gaussian approximation which should be good - // enough for our purpose. + // Tail Cancellation by Deconvolution for PASA v4 TRF // - // The parameters for the response function - const Float_t aa = 0.8872; - const Float_t bb = -0.00573; - const Float_t cc = 0.454; - const Float_t cc2 = cc*cc; + Double_t rates[2]; + Double_t coefficients[2]; - Float_t pr = aa * (bb + TMath::Exp(-x*x / (2. * cc2))); + // initialize (coefficient = alpha, rates = lambda) + + Double_t R1 = 1.0; + Double_t R2 = 1.0; + Double_t C1 = 0.5; + Double_t C2 = 0.5; + + if (nexp == 1) { // 1 Exponentials + R1 = 1.156; + R2 = 0.130; + C1 = 0.066; + C2 = 0.000; + } + if (nexp == 2) { // 2 Exponentials + R1 = 1.156; + R2 = 0.130; + C1 = 0.114; + C2 = 0.624; + } - return (pr); + 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; + + /* attention: computation order is important */ + correction=0.0; + + for ( k=0; k