X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDclusterizerV1.cxx;h=3e19c295ccbedba0e0c482886c0d52527cb0a048;hb=9679d5dbb99896e8003cf8573df706e40e172176;hp=e62926bdfd3f675cfe175ad6712657979c9f5f66;hpb=5443e65ea0e4ba7bbf3dd08d9893a31c9dd8a539;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDclusterizerV1.cxx b/TRD/AliTRDclusterizerV1.cxx index e62926bdfd3..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,69 +14,7 @@ * provided "as is" without express or implied warranty. * **************************************************************************/ -/* -$Log$ -Revision 1.16 2002/03/25 20:01:30 cblume -Introduce parameter class - -Revision 1.15 2001/11/14 12:09:11 cblume -Use correct name for digitizer - -Revision 1.14 2001/11/14 10:50:45 cblume -Changes in digits IO. Add merging of summable digits - -Revision 1.13 2001/05/28 17:07:58 hristov -Last minute changes; ExB correction in AliTRDclusterizerV1; taking into account of material in G10 TEC frames and material between TEC planes (C.Blume,S.Sedykh) - -Revision 1.12 2001/05/21 17:42:58 hristov -Constant casted to avoid the ambiguity - -Revision 1.11 2001/05/21 16:45:47 hristov -Last minute changes (C.Blume) - -Revision 1.10 2001/05/07 08:06:44 cblume -Speedup of the code. Create only AliTRDcluster - -Revision 1.9 2000/11/01 14:53:20 cblume -Merge with TRD-develop - -Revision 1.1.4.5 2000/10/15 23:40:01 cblume -Remove AliTRDconst - -Revision 1.1.4.4 2000/10/06 16:49:46 cblume -Made Getters const - -Revision 1.1.4.3 2000/10/04 16:34:58 cblume -Replace include files by forward declarations - -Revision 1.1.4.2 2000/09/22 14:49:49 cblume -Adapted to tracking code - -Revision 1.8 2000/10/02 21:28:19 fca -Removal of useless dependecies via forward declarations - -Revision 1.7 2000/06/27 13:08:50 cblume -Changed to Copy(TObject &A) to appease the HP-compiler - -Revision 1.6 2000/06/09 11:10:07 cblume -Compiler warnings and coding conventions, next round - -Revision 1.5 2000/06/08 18:32:58 cblume -Make code compliant to coding conventions - -Revision 1.4 2000/06/07 16:27:01 cblume -Try to remove compiler warnings on Sun and HP - -Revision 1.3 2000/05/08 16:17:27 cblume -Merge TRD-develop - -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$ */ /////////////////////////////////////////////////////////////////////////////// // // @@ -89,16 +28,23 @@ Add new TRD classes #include #include "AliRun.h" +#include "AliRunLoader.h" +#include "AliLoader.h" +#include "AliRawReader.h" -#include "AliTRD.h" #include "AliTRDclusterizerV1.h" #include "AliTRDmatrix.h" #include "AliTRDgeometry.h" -#include "AliTRDdigitizer.h" #include "AliTRDdataArrayF.h" #include "AliTRDdataArrayI.h" #include "AliTRDdigitsManager.h" -#include "AliTRDparameter.h" +#include "AliTRDpadPlane.h" +#include "AliTRDrawData.h" +#include "AliTRDcalibDB.h" +#include "AliTRDSimParam.h" +#include "AliTRDRecParam.h" +#include "AliTRDCommonParam.h" +#include "AliTRDcluster.h" ClassImp(AliTRDclusterizerV1) @@ -128,6 +74,7 @@ AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title //_____________________________________________________________________________ AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c) +:AliTRDclusterizer(c) { // // AliTRDclusterizerV1 copy constructor @@ -164,7 +111,7 @@ AliTRDclusterizerV1 &AliTRDclusterizerV1::operator=(const AliTRDclusterizerV1 &c } //_____________________________________________________________________________ -void AliTRDclusterizerV1::Copy(TObject &c) +void AliTRDclusterizerV1::Copy(TObject &c) const { // // Copy function @@ -183,16 +130,32 @@ Bool_t AliTRDclusterizerV1::ReadDigits() // Reads the digits arrays from the input aliroot file // - if (!fInputFile) { + if (!fRunLoader) { printf(" "); printf("No input file open\n"); return kFALSE; } - - fDigitsManager->Open(fInputFile->GetName()); + 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::ReadDigits(AliRawReader* rawReader) +{ + // + // Reads the digits arrays from the ddl file + // + + AliTRDrawData raw; + raw.SetDebug(1); + + fDigitsManager = raw.Raw2Digits(rawReader); + + return kTRUE; } @@ -205,44 +168,66 @@ Bool_t AliTRDclusterizerV1::MakeClusters() Int_t row, col, time; + /* if (fTRD->IsVersion() != 1) { printf(" "); printf("TRD must be version 1 (slow simulator).\n"); return kFALSE; } + */ // Get the geometry - AliTRDgeometry *geo = fTRD->GetGeometry(); - - // Create a default parameter class if none is defined - if (!fPar) { - fPar = new AliTRDparameter("TRDparameter","Standard TRD parameter"); + AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader); + AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); + if (!calibration) + { printf(" "); - printf("Create the default parameter object.\n"); + 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(); - Float_t timeBinSize = fPar->GetTimeBinSize(); - // Half of ampl.region - const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; - - Float_t omegaTau = fPar->GetOmegaTau(); if (fVerbose > 0) { - printf(" "); - printf("OmegaTau = %f \n",omegaTau); + //printf(" "); + //printf("OmegaTau = %f \n",omegaTau); printf(" "); printf("Start creating clusters.\n"); } - AliTRDdataArrayI *digits; + AliTRDdataArrayI *digitsIn; AliTRDdataArrayI *track0; AliTRDdataArrayI *track1; AliTRDdataArrayI *track2; // Threshold value for the maximum - Int_t maxThresh = fPar->GetClusMaxThresh(); + Float_t maxThresh = recParam->GetClusMaxThresh(); // Threshold value for the digit signal - Int_t sigThresh = fPar->GetClusSigThresh(); - + Float_t sigThresh = recParam->GetClusSigThresh(); // Iteration limit for unfolding procedure const Float_t kEpsilon = 0.01; @@ -250,52 +235,40 @@ Bool_t AliTRDclusterizerV1::MakeClusters() const Int_t kNsig = 5; const Int_t kNtrack = 3 * kNclus; - Int_t iType = 0; - Int_t iUnfold = 0; + 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]; - Float_t ratioLeft = 1.0; - Float_t ratioRight = 1.0; + 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(); - Float_t padSignal[kNsig]; - Float_t clusterSignal[kNclus]; - Float_t clusterPads[kNclus]; - Int_t clusterDigit[kNclus]; - Int_t clusterTracks[kNtrack]; + Int_t nTimeTotal = calibration->GetNumberOfTimeBins(); - Int_t chamBeg = 0; - Int_t chamEnd = AliTRDgeometry::Ncham(); - if (fTRD->GetSensChamber() >= 0) { - chamBeg = fTRD->GetSensChamber(); - chamEnd = chamBeg + 1; - } - Int_t planBeg = 0; - Int_t planEnd = AliTRDgeometry::Nplan(); - if (fTRD->GetSensPlane() >= 0) { - planBeg = fTRD->GetSensPlane(); - planEnd = planBeg + 1; + if (fVerbose > 0) { + printf(" "); + printf("Number of Time Bins = %d.\n",nTimeTotal); } - Int_t sectBeg = 0; - Int_t sectEnd = AliTRDgeometry::Nsect(); // 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++) { - if (fTRD->GetSensSector() >= 0) { - Int_t sens1 = fTRD->GetSensSector(); - Int_t sens2 = sens1 + fTRD->GetSensSectorRange(); - sens2 -= ((Int_t) (sens2 / AliTRDgeometry::Nsect())) - * AliTRDgeometry::Nsect(); - if (sens1 < sens2) { - if ((isect < sens1) || (isect >= sens2)) continue; - } - else { - if ((isect < sens1) && (isect >= sens2)) continue; - } - } + Int_t idet = geo->GetDetector(iplan,icham,isect); - Int_t idet = geo->GetDetector(iplan,icham,isect); + Int_t nRowMax = commonParam->GetRowMax(iplan,icham,isect); + Int_t nColMax = commonParam->GetColMax(iplan); Int_t nClusters = 0; Int_t nClusters2pad = 0; @@ -310,19 +283,15 @@ Bool_t AliTRDclusterizerV1::MakeClusters() ,icham,iplan,isect); } - Int_t nRowMax = fPar->GetRowMax(iplan,icham,isect); - Int_t nColMax = fPar->GetColMax(iplan); - Int_t nTimeBefore = fPar->GetTimeBefore(); - Int_t nTimeTotal = fPar->GetTimeTotal(); - - Float_t row0 = fPar->GetRow0(iplan,icham,isect); - Float_t col0 = fPar->GetCol0(iplan); - Float_t rowSize = fPar->GetRowPadSize(iplan,icham,isect); - Float_t colSize = fPar->GetColPadSize(iplan); + AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham); // Get the digits - digits = fDigitsManager->GetDigits(idet); - digits->Expand(); + 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); @@ -332,24 +301,32 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // Loop through the chamber and find the maxima for ( row = 0; row < nRowMax; row++) { - for ( col = 2; col < nColMax; col++) { + for ( col = 2; col < nColMax; col++) { + //for ( col = 4; col < nColMax-2; col++) { for (time = 0; time < nTimeTotal; time++) { - Int_t signalL = TMath::Abs(digits->GetDataUnchecked(row,col ,time)); - Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); - Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,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 (((signalL >= sigThresh) && - (signalL < signalM)) || - ((signalR >= sigThresh) && - (signalR < signalM))) { + 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 - digits->SetDataUnchecked(row,col-1,time,-signalM); + digitsOut->SetDataUnchecked(row,col-1,time,-signalM); } } - } } } @@ -360,15 +337,14 @@ Bool_t AliTRDclusterizerV1::MakeClusters() for ( col = 1; col < nColMax-1; col++) { // Maximum found ? - if (digits->GetDataUnchecked(row,col,time) < 0) { + 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(digits->GetDataUnchecked(row + clusterSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row ,iPadCol ,time)); - clusterDigit[iPad] = digits->GetIndexUnchecked(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; @@ -377,14 +353,14 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // Count the number of pads in the cluster Int_t nPadCount = 0; Int_t ii = 0; - while (TMath::Abs(digits->GetDataUnchecked(row,col-ii ,time)) + while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) { nPadCount++; ii++; if (col-ii < 0) break; } ii = 0; - while (TMath::Abs(digits->GetDataUnchecked(row,col+ii+1,time)) + while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) >= sigThresh) { nPadCount++; ii++; @@ -415,22 +391,19 @@ Bool_t AliTRDclusterizerV1::MakeClusters() break; }; - // Don't analyze large clusters - //if (iType == 4) continue; - - // Look for 5 pad cluster with minimum in the middle + // Look for 5 pad cluster with minimum in the middle Bool_t fivePadCluster = kFALSE; if (col < nColMax-3) { - if (digits->GetDataUnchecked(row,col+2,time) < 0) { + if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) { fivePadCluster = kTRUE; } if ((fivePadCluster) && (col < nColMax-5)) { - if (digits->GetDataUnchecked(row,col+4,time) >= sigThresh) { + if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) { fivePadCluster = kFALSE; } } if ((fivePadCluster) && (col > 1)) { - if (digits->GetDataUnchecked(row,col-2,time) >= sigThresh) { + if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) { fivePadCluster = kFALSE; } } @@ -441,135 +414,115 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // of the cluster which remains from a previous unfolding if (iUnfold) { clusterSignal[0] *= ratioLeft; - iType = 3; + iType = 5; iUnfold = 0; } // Unfold the 5 pad cluster if (fivePadCluster) { for (iPad = 0; iPad < kNsig; iPad++) { - padSignal[iPad] = TMath::Abs(digits->GetDataUnchecked(row - ,col-1+iPad - ,time)); + 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 = 3; + 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 - nTimeBefore + 0.5; - - if (fPar->LUTOn()) { + clusterPads[2] = time + 0.5; + + if (recParam->LUTOn()) { // Calculate the position of the cluster by using the // lookup table method - clusterPads[1] = col + 0.5 - + fPar->LUTposition(iplan,clusterSignal[0] - ,clusterSignal[1] - ,clusterSignal[2]); - + 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 - clusterPads[1] = col + 0.5 - + (clusterSignal[2] - clusterSignal[0]) - / clusterCharge; + 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); } - Float_t q0 = clusterSignal[0]; - Float_t q1 = clusterSignal[1]; - Float_t q2 = clusterSignal[2]; - Float_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) / - (clusterCharge*clusterCharge); - - // Correct for ExB displacement - if (fPar->ExBOn()) { - Int_t local_time_bin = (Int_t) clusterPads[2]; - Float_t driftLength = local_time_bin * timeBinSize + kAmWidth; - Float_t colSize = fPar->GetColPadSize(iplan); - Float_t deltaY = omegaTau*driftLength/colSize; - clusterPads[1] = clusterPads[1] - deltaY; - } - - if (fVerbose > 1) { - printf("-----------------------------------------------------------\n"); - printf("Create cluster no. %d\n",nClusters); - printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0] - ,clusterPads[1] - ,clusterPads[2]); - printf("Indices: %d, %d, %d\n",clusterDigit[0] - ,clusterDigit[1] - ,clusterDigit[2]); - printf("Total charge = %f\n",clusterCharge); - printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0] - ,clusterTracks[1] - ,clusterTracks[2]); - printf(" pad1 %d, %d, %d\n",clusterTracks[3] - ,clusterTracks[4] - ,clusterTracks[5]); - printf(" pad2 %d, %d, %d\n",clusterTracks[6] - ,clusterTracks[7] - ,clusterTracks[8]); - printf("Type = %d, Number of pads = %d\n",iType,nPadCount); - } - - // Calculate the position and the error - Float_t clusterPos[3]; - clusterPos[0] = clusterPads[1] * colSize + col0; - clusterPos[1] = clusterPads[0] * rowSize + row0; - clusterPos[2] = clusterPads[2]; - Float_t clusterSig[2]; + 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.; - + clusterSig[1] = rowSize * rowSize / 12.; + + // Add the cluster to the output array - fTRD->AddCluster(clusterPos - ,idet - ,clusterCharge - ,clusterTracks - ,clusterSig - ,iType); - + 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); } } } - } + } + + delete digitsOut; - // Compress the arrays - digits->Compress(1,0); + // Compress the arrays track0->Compress(1,0); - track1->Compress(1,0); + track1->Compress(1,0); track2->Compress(1,0); // Write the cluster and reset the array WriteClusters(idet); - fTRD->ResetRecPoints(); - - if (fVerbose > 0) { - printf(" "); - printf("Found %d clusters in total.\n" - ,nClusters); - printf(" 2pad: %d\n",nClusters2pad); - printf(" 3pad: %d\n",nClusters3pad); - printf(" 4pad: %d\n",nClusters4pad); - printf(" 5pad: %d\n",nClusters5pad); - printf(" Large: %d\n",nClustersLarge); - } - + ResetRecPoints(); } } } @@ -579,12 +532,25 @@ Bool_t AliTRDclusterizerV1::MakeClusters() printf("Done.\n"); } + //delete digitsIn; + return kTRUE; } //_____________________________________________________________________________ -Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, 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. @@ -593,15 +559,23 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, Float_t* padSignal // The resulting ratio is then returned to the calling method. // - Int_t irc = 0; - 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 - Float_t newSignal[3] = {0}; + 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 the iteration while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { @@ -610,23 +584,23 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, Float_t* padSignal 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]); + 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]); // Set cluster charge ratio - irc = fPar->PadResponse(1.0,maxLeft ,plane,newSignal); - Float_t ampLeft = padSignal[1] / newSignal[1]; - irc = fPar->PadResponse(1.0,maxRight,plane,newSignal); - Float_t ampRight = padSignal[3] / newSignal[1]; + 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]; // Apply pad response to parameters - irc = fPar->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal ); - irc = fPar->PadResponse(ampRight,maxRight,plane,newRightSignal); + irc = calibration->PadResponse(ampLeft ,maxLeft ,plane,newLeftSignal ); + irc = calibration->PadResponse(ampRight,maxRight,plane,newRightSignal); // Calculate new overlapping ratio - ratio = TMath::Min((Float_t)1.0,newLeftSignal[2] / + ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] / (newLeftSignal[2] + newRightSignal[0])); } @@ -635,3 +609,141 @@ Float_t AliTRDclusterizerV1::Unfold(Float_t eps, Int_t plane, Float_t* padSignal } +//_____________________________________________________________________________ +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) +{ + // + // Tail Cancellation by Deconvolution for PASA v4 TRF + // + + Double_t rates[2]; + Double_t coefficients[2]; + + // 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; + } + + 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