X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TRD%2FAliTRDclusterizerV1.cxx;h=4b7f245f95f0fb936b9f03bd35bc4507a408b115;hb=9bce4ac18e166984365f8b31b0c44d05d8de6fb3;hp=3e19c295ccbedba0e0c482886c0d52527cb0a048;hpb=4ab6879608b6572fe1e06a2c438408d50a93c212;p=u%2Fmrichter%2FAliRoot.git diff --git a/TRD/AliTRDclusterizerV1.cxx b/TRD/AliTRDclusterizerV1.cxx index 3e19c295ccb..4b7f245f95f 100644 --- a/TRD/AliTRDclusterizerV1.cxx +++ b/TRD/AliTRDclusterizerV1.cxx @@ -18,7 +18,7 @@ /////////////////////////////////////////////////////////////////////////////// // // -// TRD cluster finder for the slow simulator. +// TRD cluster finder // // // /////////////////////////////////////////////////////////////////////////////// @@ -27,13 +27,13 @@ #include #include -#include "AliRun.h" #include "AliRunLoader.h" #include "AliLoader.h" #include "AliRawReader.h" +#include "AliLog.h" +#include "AliAlignObj.h" #include "AliTRDclusterizerV1.h" -#include "AliTRDmatrix.h" #include "AliTRDgeometry.h" #include "AliTRDdataArrayF.h" #include "AliTRDdataArrayI.h" @@ -43,45 +43,46 @@ #include "AliTRDcalibDB.h" #include "AliTRDSimParam.h" #include "AliTRDRecParam.h" -#include "AliTRDCommonParam.h" #include "AliTRDcluster.h" +#include "Cal/AliTRDCalROC.h" +#include "Cal/AliTRDCalDet.h" + ClassImp(AliTRDclusterizerV1) //_____________________________________________________________________________ -AliTRDclusterizerV1::AliTRDclusterizerV1():AliTRDclusterizer() +AliTRDclusterizerV1::AliTRDclusterizerV1() + :AliTRDclusterizer() + ,fDigitsManager(NULL) { // // AliTRDclusterizerV1 default constructor // - fDigitsManager = 0; - } //_____________________________________________________________________________ -AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t* name, const Text_t* title) - :AliTRDclusterizer(name,title) +AliTRDclusterizerV1::AliTRDclusterizerV1(const Text_t *name, const Text_t *title) + :AliTRDclusterizer(name,title) + ,fDigitsManager(new AliTRDdigitsManager()) { // - // AliTRDclusterizerV1 default constructor + // AliTRDclusterizerV1 constructor // - fDigitsManager = new AliTRDdigitsManager(); fDigitsManager->CreateArrays(); } //_____________________________________________________________________________ AliTRDclusterizerV1::AliTRDclusterizerV1(const AliTRDclusterizerV1 &c) -:AliTRDclusterizer(c) + :AliTRDclusterizer(c) + ,fDigitsManager(NULL) { // // AliTRDclusterizerV1 copy constructor // - ((AliTRDclusterizerV1 &) c).Copy(*this); - } //_____________________________________________________________________________ @@ -131,12 +132,14 @@ Bool_t AliTRDclusterizerV1::ReadDigits() // if (!fRunLoader) { - printf(" "); - printf("No input file open\n"); + AliError("No run loader available"); return kFALSE; } + AliLoader* loader = fRunLoader->GetLoader("TRDLoader"); - if (!loader->TreeD()) loader->LoadDigits(); + if (!loader->TreeD()) { + loader->LoadDigits(); + } // Read in the digit arrays return (fDigitsManager->ReadDigits(loader->TreeD())); @@ -144,15 +147,25 @@ Bool_t AliTRDclusterizerV1::ReadDigits() } //_____________________________________________________________________________ -Bool_t AliTRDclusterizerV1::ReadDigits(AliRawReader* rawReader) +Bool_t AliTRDclusterizerV1::ReadDigits(TTree *digitsTree) +{ + // + // Reads the digits arrays from the input tree + // + + // Read in the digit arrays + return (fDigitsManager->ReadDigits(digitsTree)); + +} + +//_____________________________________________________________________________ +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; @@ -166,170 +179,160 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // Generates the cluster. // - Int_t row, col, time; + Int_t row = 0; + Int_t col = 0; + Int_t time = 0; + Int_t icham = 0; + Int_t iplan = 0; + Int_t isect = 0; + Int_t iPad = 0; + + AliTRDdataArrayI *digitsIn; + AliTRDdataArrayI *tracksIn; - /* - if (fTRD->IsVersion() != 1) { - printf(" "); - printf("TRD must be version 1 (slow simulator).\n"); - return kFALSE; - } - */ - - // Get the geometry - AliTRDgeometry *geo = AliTRDgeometry::GetGeometry(fRunLoader); - AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); - if (!calibration) - { - printf(" "); - printf("ERROR getting instance of AliTRDcalibDB"); + AliTRDgeometry geo; + + AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); + if (!calibration) { + AliFatal("No AliTRDcalibDB instance available\n"); return kFALSE; } - AliTRDSimParam* simParam = AliTRDSimParam::Instance(); - if (!simParam) - { - printf(" "); - printf("ERROR getting instance of AliTRDSimParam"); + AliTRDSimParam *simParam = AliTRDSimParam::Instance(); + if (!simParam) { + AliError("No AliTRDSimParam instance available\n"); return kFALSE; } - AliTRDRecParam* recParam = AliTRDRecParam::Instance(); - if (!recParam) - { - printf(" "); - printf("ERROR getting instance of AliTRDRecParam"); + AliTRDRecParam *recParam = AliTRDRecParam::Instance(); + if (!recParam) { + AliError("No AliTRDRecParam instance available\n"); 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; + // ADC thresholds + // Float_t ADCthreshold = simParam->GetADCthreshold(); + Float_t ADCthreshold = 0; // Threshold value for the maximum - Float_t maxThresh = recParam->GetClusMaxThresh(); + Float_t maxThresh = recParam->GetClusMaxThresh(); // Threshold value for the digit signal - Float_t sigThresh = recParam->GetClusSigThresh(); + Float_t sigThresh = recParam->GetClusSigThresh(); + + // Detector wise calibration object for t0 + const AliTRDCalDet *calT0Det = calibration->GetT0Det(); + // Detector wise calibration object for the gain factors + const AliTRDCalDet *calGainFactorDet = calibration->GetGainFactorDet(); + // Iteration limit for unfolding procedure const Float_t kEpsilon = 0.01; - const Int_t kNclus = 3; const Int_t kNsig = 5; - const Int_t kNtrack = 3 * kNclus; + const Int_t kNdict = AliTRDdigitsManager::kNDict; + const Int_t kNtrack = kNdict * kNclus; - Int_t iType = 0; Int_t iUnfold = 0; Double_t ratioLeft = 1.0; Double_t ratioRight = 1.0; - // + Int_t iClusterROC = 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 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(); - Int_t nTimeTotal = calibration->GetNumberOfTimeBins(); - - if (fVerbose > 0) { - printf(" "); - printf("Number of Time Bins = %d.\n",nTimeTotal); - } + AliDebug(1,Form("Number of Time Bins = %d.\n",nTimeTotal)); // 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); + for (icham = chamBeg; icham < chamEnd; icham++) { + for (iplan = planBeg; iplan < planEnd; iplan++) { + for (isect = sectBeg; isect < sectEnd; isect++) { - Int_t nRowMax = commonParam->GetRowMax(iplan,icham,isect); - Int_t nColMax = commonParam->GetColMax(iplan); - - Int_t nClusters = 0; - Int_t nClusters2pad = 0; - Int_t nClusters3pad = 0; - Int_t nClusters4pad = 0; - Int_t nClusters5pad = 0; - Int_t nClustersLarge = 0; - - if (fVerbose > 0) { - printf(" "); - printf("Analyzing chamber %d, plane %d, sector %d.\n" - ,icham,iplan,isect); - } - - AliTRDpadPlane *padPlane = commonParam->GetPadPlane(iplan,icham); + Int_t idet = geo.GetDetector(iplan,icham,isect); + Int_t ilayer = AliGeomManager::kTRD1 + iplan; + Int_t imodule = icham + chamEnd * isect; + UShort_t volid = AliGeomManager::LayerToVolUID(ilayer,imodule); // Get the digits digitsIn = fDigitsManager->GetDigits(idet); + // This is to take care of switched off super modules + if (digitsIn->GetNtime() == 0) { + continue; + } digitsIn->Expand(); - AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow(), digitsIn->GetNcol(), digitsIn->GetNtime()); + AliTRDdataArrayI *tracksTmp = fDigitsManager->GetDictionary(idet,0); + tracksTmp->Expand(); + + Int_t nRowMax = geo.GetRowMax(iplan,icham,isect); + Int_t nColMax = geo.GetColMax(iplan); + + AliTRDpadPlane *padPlane = geo.GetPadPlane(iplan,icham); - Transform(digitsIn, digitsOut, idet, nRowMax, nColMax, nTimeTotal, ADCthreshold); + // Calibration object with pad wise values for t0 + AliTRDCalROC *calT0ROC = calibration->GetT0ROC(idet); + // Calibration object with pad wise values for the gain factors + AliTRDCalROC *calGainFactorROC = calibration->GetGainFactorROC(idet); + // Calibration value for chamber wise t0 + Float_t calT0DetValue = calT0Det->GetValue(idet); + // Calibration value for chamber wise gain factor + Float_t calGainFactorDetValue = calGainFactorDet->GetValue(idet); - track0 = fDigitsManager->GetDictionary(idet,0); - track0->Expand(); - track1 = fDigitsManager->GetDictionary(idet,1); - track1->Expand(); - track2 = fDigitsManager->GetDictionary(idet,2); - track2->Expand(); + Int_t nClusters = 0; + + // Apply the gain and the tail cancelation via digital filter + AliTRDdataArrayF *digitsOut = new AliTRDdataArrayF(digitsIn->GetNrow() + ,digitsIn->GetNcol() + ,digitsIn->GetNtime()); + Transform(digitsIn + ,digitsOut + ,nRowMax,nColMax,nTimeTotal + ,ADCthreshold + ,calGainFactorROC + ,calGainFactorDetValue); + + // Input digits are not needed any more + digitsIn->Compress(1,0); // 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); + + Float_t signalL = TMath::Abs(digitsOut->GetDataUnchecked(row,col ,time)); + Float_t signalR = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)); + + if ((TMath::Abs(signalL) <= signalM) && + (TMath::Abs(signalR) < signalM)) { + if ((TMath::Abs(signalL) >= sigThresh) || + (TMath::Abs(signalR) >= sigThresh)) { + // Maximum found, mark the position by a negative signal + digitsOut->SetDataUnchecked(row,col-1,time,-signalM); + } } + } - } - } - } + + } + } + } + tracksTmp->Compress(1,0); + + // The index to the first cluster of a given ROC + Int_t firstClusterROC = -1; + // The number of cluster in a given ROC + Int_t nClusterROC = 0; // Now check the maxima and calculate the cluster position for ( row = 0; row < nRowMax ; row++) { @@ -337,72 +340,45 @@ Bool_t AliTRDclusterizerV1::MakeClusters() for ( col = 1; col < nColMax-1; col++) { // Maximum found ? - if (digitsOut->GetDataUnchecked(row,col,time) < 0) { + if (digitsOut->GetDataUnchecked(row,col,time) < 0.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; + clusterSignal[iPad] = + TMath::Abs(digitsOut->GetDataUnchecked(row,iPadCol,time)); } // 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) { + Int_t ii; + // Look to the left + ii = 0; + while (TMath::Abs(digitsOut->GetDataUnchecked(row,col-ii ,time)) >= sigThresh) { nPadCount++; ii++; if (col-ii < 0) break; } + // Look to the right ii = 0; - while (TMath::Abs(digitsOut->GetDataUnchecked(row,col+ii+1,time)) - >= sigThresh) { + 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 + + // Look for 5 pad cluster with minimum in the middle Bool_t fivePadCluster = kFALSE; - if (col < nColMax-3) { + if (col < (nColMax - 3)) { if (digitsOut->GetDataUnchecked(row,col+2,time) < 0) { fivePadCluster = kTRUE; } - if ((fivePadCluster) && (col < nColMax-5)) { + if ((fivePadCluster) && (col < (nColMax - 5))) { if (digitsOut->GetDataUnchecked(row,col+4,time) >= sigThresh) { fivePadCluster = kFALSE; } } - if ((fivePadCluster) && (col > 1)) { + if ((fivePadCluster) && (col > 1)) { if (digitsOut->GetDataUnchecked(row,col-2,time) >= sigThresh) { fivePadCluster = kFALSE; } @@ -414,7 +390,6 @@ Bool_t AliTRDclusterizerV1::MakeClusters() // of the cluster which remains from a previous unfolding if (iUnfold) { clusterSignal[0] *= ratioLeft; - iType = 5; iUnfold = 0; } @@ -422,15 +397,14 @@ Bool_t AliTRDclusterizerV1::MakeClusters() if (fivePadCluster) { for (iPad = 0; iPad < kNsig; iPad++) { padSignal[iPad] = TMath::Abs(digitsOut->GetDataUnchecked(row - ,col-1+iPad - ,time)); + ,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; } @@ -439,100 +413,165 @@ Bool_t AliTRDclusterizerV1::MakeClusters() + clusterSignal[2]; // The position of the cluster - clusterPads[0] = row + 0.5; + 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]); + ,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 ,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)) < padSignal[1])) { padSignal[0] = TMath::Abs(digitsOut->GetDataUnchecked(row,col-2,time)); } - if (colGetDataUnchecked(row,col+2,time)GetDataUnchecked(row,col+2,time)) < padSignal[3])) { padSignal[4] = TMath::Abs(digitsOut->GetDataUnchecked(row,col+2,time)); } - clusterPads[1] = GetCOG(padSignal); - + 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); + Double_t clusterSigmaY2 = (q1 * (q0 + q2) + 4.0 * q0 * q2) + / (clusterCharge*clusterCharge); - - + // // Calculate the position and the error + // + + // Correct for t0 (sum of chamber and pad wise values !!!) + Float_t calT0ROCValue = calT0ROC->GetValue(col,row); + Char_t clusterTimeBin = ((Char_t) TMath::Nint(time - (calT0DetValue + calT0ROCValue))); + Double_t colSize = padPlane->GetColSize(col); + Double_t rowSize = padPlane->GetRowSize(row); + + Float_t clusterPos[3]; + clusterPos[0] = padPlane->GetColPos(col) - (clusterPads[1] + 0.5) * colSize; + clusterPos[1] = padPlane->GetRowPos(row) - 0.5 * rowSize; + clusterPos[2] = CalcXposFromTimebin(clusterPads[2],idet,col,row); + Float_t clusterSig[2]; + clusterSig[0] = (clusterSigmaY2 + 1.0/12.0) * colSize*colSize; + clusterSig[1] = rowSize * rowSize / 12.0; - // 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 - 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; + // Store the amplitudes of the pads in the cluster for later analysis + 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); - } - } - } - } + + // Add the cluster to the output array + // The track indices will be stored later + AliTRDcluster *cluster = new AliTRDcluster(idet + ,clusterCharge + ,clusterPos + ,clusterSig + ,0x0 + ,((Char_t) nPadCount) + ,signals + ,((UChar_t) col) + ,clusterTimeBin + ,clusterPads[1] + ,volid); + // Temporarily store the row, column and time bin of the center pad + // Used to later on assign the track indices + cluster->SetLabel( row,0); + cluster->SetLabel( col,1); + cluster->SetLabel(time,2); + RecPoints()->Add(cluster); + + // Store the index of the first cluster in the current ROC + if (firstClusterROC < 0) { + firstClusterROC = RecPoints()->GetEntriesFast() - 1; + } + // Count the number of cluster in the current ROC + nClusterROC++; + + } // if: Maximum found ? + + } // loop: pad columns + } // loop: time bins + } // loop: pad rows delete digitsOut; - // Compress the arrays - track0->Compress(1,0); - track1->Compress(1,0); - track2->Compress(1,0); + // + // Add the track indices to the found clusters + // + + // Temporary array to collect the track indices + Int_t *idxTracks = new Int_t[kNtrack*nClusterROC]; + + // Loop through the dictionary arrays one-by-one + // to keep memory consumption low + for (Int_t iDict = 0; iDict < kNdict; iDict++) { + + tracksIn = fDigitsManager->GetDictionary(idet,iDict); + tracksIn->Expand(); + + // Loop though the clusters found in this ROC + for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) { + + AliTRDcluster *cluster = (AliTRDcluster *) + RecPoints()->UncheckedAt(firstClusterROC+iClusterROC); + row = cluster->GetLabel(0); + col = cluster->GetLabel(1); + time = cluster->GetLabel(2); + + for (iPad = 0; iPad < kNclus; iPad++) { + Int_t iPadCol = col - 1 + iPad; + Int_t index = tracksIn->GetDataUnchecked(row,iPadCol,time) - 1; + idxTracks[3*iPad+iDict + iClusterROC*kNtrack] = index; + } + + } + + // Compress the arrays + tracksIn->Compress(1,0); + + } + + // Copy the track indices into the cluster + // Loop though the clusters found in this ROC + for (iClusterROC = 0; iClusterROC < nClusterROC; iClusterROC++) { + + AliTRDcluster *cluster = (AliTRDcluster *) + RecPoints()->UncheckedAt(firstClusterROC+iClusterROC); + cluster->SetLabel(-9999,0); + cluster->SetLabel(-9999,1); + cluster->SetLabel(-9999,2); + + cluster->AddTrackIndex(&idxTracks[iClusterROC*kNtrack]); + + } + + delete [] idxTracks; // Write the cluster and reset the array WriteClusters(idet); ResetRecPoints(); - } - } - } - - if (fVerbose > 0) { - printf(" "); - printf("Done.\n"); - } - //delete digitsIn; + } // loop: Sectors + } // loop: Planes + } // loop: Chambers return kTRUE; @@ -542,15 +581,25 @@ Bool_t AliTRDclusterizerV1::MakeClusters() 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; + // 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) +Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t *padSignal) { // // Method to unfold neighbouring maxima. @@ -559,23 +608,21 @@ Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSig // The resulting ratio is then returned to the calling method. // - AliTRDcalibDB* calibration = AliTRDcalibDB::Instance(); - if (!calibration) - { - printf(" "); - printf("ERROR getting instance of AliTRDcalibDB"); + AliTRDcalibDB *calibration = AliTRDcalibDB::Instance(); + if (!calibration) { + AliError("No AliTRDcalibDB instance available\n"); return kFALSE; } Int_t irc = 0; - Int_t itStep = 0; // Count iteration steps + Int_t itStep = 0; // Count iteration steps - Double_t ratio = 0.5; // Start value for ratio - Double_t prevRatio = 0; // Store previous ratio + Double_t ratio = 0.5; // Start value for ratio + Double_t prevRatio = 0.0; // Store previous ratio - 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}; + Double_t newLeftSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store left cluster signal + Double_t newRightSignal[3] = { 0.0, 0.0, 0.0 }; // Array to store right cluster signal + Double_t newSignal[3] = { 0.0, 0.0, 0.0 }; // Start the iteration while ((TMath::Abs(prevRatio - ratio) > eps) && (itStep < 10)) { @@ -587,7 +634,7 @@ Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSig 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]); + / ((1.0 - ratio)*padSignal[2] + padSignal[3] + padSignal[4]); // Set cluster charge ratio irc = calibration->PadResponse(1.0,maxLeft ,plane,newSignal); @@ -601,7 +648,7 @@ Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSig // Calculate new overlapping ratio ratio = TMath::Min((Double_t)1.0,newLeftSignal[2] / - (newLeftSignal[2] + newRightSignal[0])); + (newLeftSignal[2] + newRightSignal[0])); } @@ -610,74 +657,64 @@ Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSig } //_____________________________________________________________________________ -void AliTRDclusterizerV1::Transform(AliTRDdataArrayI* digitsIn, - AliTRDdataArrayF* digitsOut, - Int_t idet, Int_t nRowMax, - Int_t nColMax, Int_t nTimeTotal, - Float_t ADCthreshold) +void AliTRDclusterizerV1::Transform(AliTRDdataArrayI *digitsIn + , AliTRDdataArrayF *digitsOut + , Int_t nRowMax, Int_t nColMax, Int_t nTimeTotal + , Float_t ADCthreshold + , AliTRDCalROC *calGainFactorROC + , Float_t calGainFactorDetValue) { - // // Apply gain factor - // Apply tail cancellation: Transform digitsIn to digitsOut + // Apply tail cancelation: Transform digitsIn to digitsOut // + Int_t iRow = 0; + Int_t iCol = 0; + Int_t iTime = 0; - AliTRDRecParam* recParam = AliTRDRecParam::Instance(); - if (!recParam) - { - printf(" "); - printf("ERROR getting instance of AliTRDRecParam"); + AliTRDRecParam *recParam = AliTRDRecParam::Instance(); + if (!recParam) { + AliError("No AliTRDRecParam instance available\n"); 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 + 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 (iRow = 0; iRow < nRowMax; iRow++ ) { + for (iCol = 0; iCol < nColMax; iCol++ ) { + + Float_t calGainFactorROCValue = calGainFactorROC->GetValue(iCol,iRow); + Double_t gain = calGainFactorDetValue + * calGainFactorROCValue; + + for (iTime = 0; iTime < nTimeTotal; iTime++) { - 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 + // 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] = digitsIn->GetDataUnchecked(iRow,iCol,iTime); + inADC[iTime] /= gain; + outADC[iTime] = inADC[iTime]; - inADC[iTime] /= gain; - outADC[iTime] = inADC[iTime]; } // Apply the tail cancelation via the digital filter - if (recParam->TCOn()) - { + if (recParam->TCOn()) { DeConvExp(inADC,outADC,nTimeTotal,recParam->GetTCnexp()); } - for (Int_t iTime = 0; iTime < nTimeTotal; iTime++) { + for (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; @@ -687,20 +724,18 @@ void AliTRDclusterizerV1::Transform(AliTRDdataArrayI* digitsIn, } - //_____________________________________________________________________________ -void AliTRDclusterizerV1::DeConvExp(Double_t *source, Double_t *target, - Int_t n, Int_t nexp) +void AliTRDclusterizerV1::DeConvExp(Double_t *source, Double_t *target + , Int_t n, Int_t nexp) { // - // Tail Cancellation by Deconvolution for PASA v4 TRF + // Tail cancellation by deconvolution for PASA v4 TRF // Double_t rates[2]; Double_t coefficients[2]; - // initialize (coefficient = alpha, rates = lambda) - + // Initialization (coefficient = alpha, rates = lambda) Double_t R1 = 1.0; Double_t R2 = 1.0; Double_t C1 = 0.5; @@ -722,28 +757,34 @@ void AliTRDclusterizerV1::DeConvExp(Double_t *source, Double_t *target, coefficients[0] = C1; coefficients[1] = C2; - Double_t Dt = 0.100; + Double_t Dt = 0.1; rates[0] = TMath::Exp(-Dt/(R1)); rates[1] = TMath::Exp(-Dt/(R2)); - Int_t i, k; - Double_t reminder[2]; - Double_t correction, result; + Int_t i = 0; + Int_t k = 0; - /* attention: computation order is important */ - correction=0.0; - - for ( k=0; k