X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCclustererMI.cxx;h=769730b547da36f6877c8625d4765d6c61d5f1d3;hb=ce8cd162576831b6cef9b8b57652c77af22b740a;hp=d546627eff1dc3d2fcd6eecbc43a6461bcfefdcb;hpb=d7a1155515f6473e31d47feb1d2c09bd2b63aa11;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCclustererMI.cxx b/TPC/AliTPCclustererMI.cxx index d546627eff1..769730b547d 100644 --- a/TPC/AliTPCclustererMI.cxx +++ b/TPC/AliTPCclustererMI.cxx @@ -18,37 +18,225 @@ //------------------------------------------------------- // Implementation of the TPC clusterer // +// 1. The Input data for reconstruction - Options +// 1.a Simulated data - TTree - invoked Digits2Clusters() +// 1.b Raw data - Digits2Clusters(AliRawReader* rawReader); +// +// 2. The Output data +// 2.a TTree with clusters - if SetOutput(TTree * tree) invoked +// 2.b TObjArray - Faster option for HLT +// 2.c TClonesArray - Faster option for HLT (smaller memory consumption), activate with fBClonesArray flag +// +// 3. Reconstruction setup +// see AliTPCRecoParam for list of parameters +// The reconstruction parameterization taken from the +// AliTPCReconstructor::GetRecoParam() +// Possible to setup it in reconstruction macro AliTPCReconstructor::SetRecoParam(...) +// +// +// // Origin: Marian Ivanov //------------------------------------------------------- -#include "AliTPCReconstructor.h" -#include "AliTPCclustererMI.h" -#include "AliTPCclusterMI.h" -#include +#include "Riostream.h" +#include #include -#include "AliTPCClustersArray.h" -#include "AliTPCClustersRow.h" -#include "AliTPCRawStream.h" +#include +#include +#include +#include +#include +#include +#include + #include "AliDigits.h" +#include "AliLoader.h" +#include "AliLog.h" +#include "AliMathBase.h" +#include "AliRawEventHeaderBase.h" +#include "AliRawReader.h" +#include "AliRunLoader.h" #include "AliSimDigits.h" +#include "AliTPCCalPad.h" +#include "AliTPCCalROC.h" +#include "AliTPCClustersArray.h" +#include "AliTPCClustersRow.h" #include "AliTPCParam.h" -#include "AliRawReader.h" #include "AliTPCRawStream.h" -#include "AliRunLoader.h" -#include "AliLoader.h" -#include "Riostream.h" -#include +#include "AliTPCRawStreamV3.h" +#include "AliTPCRecoParam.h" +#include "AliTPCReconstructor.h" +#include "AliTPCcalibDB.h" +#include "AliTPCclusterInfo.h" +#include "AliTPCclusterMI.h" +#include "AliTPCTransform.h" +#include "AliTPCclustererMI.h" ClassImp(AliTPCclustererMI) -AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par) +AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam): + fBins(0), + fSigBins(0), + fNSigBins(0), + fLoop(0), + fMaxBin(0), + fMaxTime(1006), // 1000>940 so use 1000, add 3 virtual time bins before and 3 after + fMaxPad(0), + fSector(-1), + fRow(-1), + fSign(0), + fRx(0), + fPadWidth(0), + fPadLength(0), + fZWidth(0), + fPedSubtraction(kFALSE), + fEventHeader(0), + fTimeStamp(0), + fEventType(0), + fInput(0), + fOutput(0), + fOutputArray(0), + fOutputClonesArray(0), + fRowCl(0), + fRowDig(0), + fParam(0), + fNcluster(0), + fNclusters(0), + fDebugStreamer(0), + fRecoParam(0), + fBDumpSignal(kFALSE), + fBClonesArray(kFALSE), + fAllBins(NULL), + fAllSigBins(NULL), + fAllNSigBins(NULL) { + // + // COSNTRUCTOR + // param - tpc parameters for given file + // recoparam - reconstruction parameters + // fInput =0; - fOutput=0; fParam = par; + if (recoParam) { + fRecoParam = recoParam; + }else{ + //set default parameters if not specified + fRecoParam = AliTPCReconstructor::GetRecoParam(); + if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam(); + } + + if(AliTPCReconstructor::StreamLevel()>0) { + fDebugStreamer = new TTreeSRedirector("TPCsignal.root"); + } + + // Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin(); + fRowCl= new AliTPCClustersRow("AliTPCclusterMI"); + + // Non-persistent arrays + // + //alocate memory for sector - maximal case + // + AliTPCROC * roc = AliTPCROC::Instance(); + Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1); + Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1); + + fAllBins = new Float_t*[nRowsMax]; + fAllSigBins = new Int_t*[nRowsMax]; + fAllNSigBins = new Int_t[nRowsMax]; + for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { + // + Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after + fAllBins[iRow] = new Float_t[maxBin]; + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + fAllSigBins[iRow] = new Int_t[maxBin]; + fAllNSigBins[iRow]=0; + } +} +//______________________________________________________________ +AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI ¶m) + :TObject(param), + fBins(0), + fSigBins(0), + fNSigBins(0), + fLoop(0), + fMaxBin(0), + fMaxTime(0), + fMaxPad(0), + fSector(-1), + fRow(-1), + fSign(0), + fRx(0), + fPadWidth(0), + fPadLength(0), + fZWidth(0), + fPedSubtraction(kFALSE), + fEventHeader(0), + fTimeStamp(0), + fEventType(0), + fInput(0), + fOutput(0), + fOutputArray(0), + fOutputClonesArray(0), + fRowCl(0), + fRowDig(0), + fParam(0), + fNcluster(0), + fNclusters(0), + fDebugStreamer(0), + fRecoParam(0), + fBDumpSignal(kFALSE), + fBClonesArray(kFALSE), + fAllBins(NULL), + fAllSigBins(NULL), + fAllNSigBins(NULL) +{ + // + // dummy + // + fMaxBin = param.fMaxBin; } +//______________________________________________________________ +AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & param) +{ + // + // assignment operator - dummy + // + fMaxBin=param.fMaxBin; + return (*this); +} +//______________________________________________________________ +AliTPCclustererMI::~AliTPCclustererMI(){ + // + // + // + if (fDebugStreamer) delete fDebugStreamer; + if (fOutputArray){ + //fOutputArray->Delete(); + delete fOutputArray; + } + if (fOutputClonesArray){ + fOutputClonesArray->Delete(); + delete fOutputClonesArray; + } + + if (fRowCl) { + fRowCl->GetArray()->Delete(); + delete fRowCl; + } + + AliTPCROC * roc = AliTPCROC::Instance(); + Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1); + for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { + delete [] fAllBins[iRow]; + delete [] fAllSigBins[iRow]; + } + delete [] fAllBins; + delete [] fAllSigBins; + delete [] fAllNSigBins; +} + void AliTPCclustererMI::SetInput(TTree * tree) { // @@ -65,19 +253,35 @@ void AliTPCclustererMI::SetInput(TTree * tree) void AliTPCclustererMI::SetOutput(TTree * tree) { // + // Set the output tree + // If not set the ObjArray used - Option for HLT // - fOutput= tree; - AliTPCClustersRow clrow; - AliTPCClustersRow *pclrow=&clrow; - clrow.SetClass("AliTPCclusterMI"); - clrow.SetArray(1); // to make Clones array + if (!tree) return; + fOutput= tree; + AliTPCClustersRow clrow, *pclrow=&clrow; + pclrow = new AliTPCClustersRow("AliTPCclusterMI"); fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200); } +void AliTPCclustererMI::FillRow(){ + // + // fill the output container - + // 2 Options possible + // Tree + // TObjArray + // + if (fOutput) fOutput->Fill(); + if (!fOutput && !fBClonesArray){ + // + if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal()); + if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID()); + } +} + Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){ // sigma y2 = in digits - we don't know the angle - Float_t z = iz*fParam->GetZWidth(); + Float_t z = iz*fParam->GetZWidth()+fParam->GetNTBinsL1()*fParam->GetZWidth(); Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/ (fPadWidth*fPadWidth); Float_t sres = 0.25; @@ -88,9 +292,9 @@ Float_t AliTPCclustererMI::GetSigmaY2(Int_t iz){ Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){ //sigma z2 = in digits - angle estimated supposing vertex constraint - Float_t z = iz*fZWidth; + Float_t z = iz*fZWidth+fParam->GetNTBinsL1()*fParam->GetZWidth(); Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth); - Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth); + Float_t angular = fPadLength*(fParam->GetZLength(fSector)-z)/(fRx*fZWidth); angular*=angular; angular/=12.; Float_t sres = fParam->GetZSigma()/fZWidth; @@ -99,19 +303,23 @@ Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){ return res; } -void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Int_t *bins, UInt_t /*m*/, +void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Float_t *bins, UInt_t /*m*/, AliTPCclusterMI &c) { + // + // k - Make cluster at position k + // bins - 2 D array of signals mapped to 1 dimensional array - + // max - the number of time bins er one dimension + // c - refernce to cluster to be filled + // Int_t i0=k/max; //central pad Int_t j0=k%max; //central time bin // set pointers to data //Int_t dummy[5] ={0,0,0,0,0}; - Int_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2 - Int_t * resmatrix[5]; + Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2 for (Int_t di=-2;di<=2;di++){ matrix[di+2] = &bins[k+di*max]; - resmatrix[di+2] = &fResBins[k+di*max]; } //build matrix with virtual charge Float_t sigmay2= GetSigmaY2(j0); @@ -120,7 +328,7 @@ AliTPCclusterMI &c) Float_t vmatrix[5][5]; vmatrix[2][2] = matrix[2][0]; c.SetType(0); - c.SetMax(Short_t(vmatrix[2][2])); // write maximal amplitude + c.SetMax((UShort_t)(vmatrix[2][2])); // write maximal amplitude for (Int_t di =-1;di <=1;di++) for (Int_t dj =-1;dj <=1;dj++){ Float_t amp = matrix[di+2][dj]; @@ -188,8 +396,10 @@ AliTPCclusterMI &c) // if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return; - if ( (ry <1.2) && (rz<1.2) ) { - //if cluster looks like expected + if ( ((ry <1.2) && (rz<1.2)) || (!fRecoParam->GetDoUnfold())) { + // + //if cluster looks like expected or Unfolding not switched on + //standard COG is used //+1.2 deviation from expected sigma accepted // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2)); @@ -197,26 +407,20 @@ AliTPCclusterMI &c) meanj +=j0; //set cluster parameters c.SetQ(sumw); - c.SetY(meani*fPadWidth); - c.SetZ(meanj*fZWidth); + c.SetPad(meani-2.5); + c.SetTimeBin(meanj-3); c.SetSigmaY2(mi2); c.SetSigmaZ2(mj2); - AddCluster(c); - //remove cluster data from data - for (Int_t di=-2;di<=2;di++) - for (Int_t dj=-2;dj<=2;dj++){ - resmatrix[di+2][dj] -= Int_t(vmatrix[di+2][dj+2]); - if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0; - } - resmatrix[2][0] =0; + c.SetType(0); + AddCluster(c,(Float_t*)vmatrix,k); return; } // //unfolding when neccessary // - Int_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3 - Int_t dummy[7]={0,0,0,0,0,0}; + Float_t * matrix2[7]; //7x7 matrix with digits - indexing i = 0 ..6 j = -3..3 + Float_t dummy[7]={0,0,0,0,0,0}; for (Int_t di=-3;di<=3;di++){ matrix2[di+3] = &bins[k+di*max]; if ((k+di*max)<3) matrix2[di+3] = &dummy[3]; @@ -232,23 +436,21 @@ AliTPCclusterMI &c) meanj +=j0; //set cluster parameters c.SetQ(sumu); - c.SetY(meani*fPadWidth); - c.SetZ(meanj*fZWidth); + c.SetPad(meani-2.5); + c.SetTimeBin(meanj-3); c.SetSigmaY2(mi2); c.SetSigmaZ2(mj2); c.SetType(Char_t(overlap)+1); - AddCluster(c); + AddCluster(c,(Float_t*)vmatrix,k); //unfolding 2 meani-=i0; meanj-=j0; - if (gDebug>4) - printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]); } -void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, +void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, Float_t & sumu, Float_t & overlap ) { // @@ -280,7 +482,7 @@ void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5 else{ Float_t ratio =1; if ( ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))|| - sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){ + (sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 )){ Float_t xm2 = sum3i[-dk+3]; Float_t xm1 = sum3i[+3]; Float_t x1 = sum3i[2*dk+3]; @@ -371,8 +573,6 @@ void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5 } } } - if (gDebug>4) - printf("%f\n", recmatrix[2][2]); } @@ -395,50 +595,94 @@ Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, F return max; } -void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){ +void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * /*matrix*/, Int_t /*pos*/){ // - // transform cluster to the global coordinata - // add the cluster to the array // - Float_t meani = c.GetY()/fPadWidth; - Float_t meanj = c.GetZ()/fZWidth; + // Transform cluster to the rotated global coordinata + // Assign labels to the cluster + // add the cluster to the array + // for more details - See AliTPCTranform::Transform(x,i,0,1) + Float_t meani = c.GetPad(); + Float_t meanj = c.GetTimeBin(); - Int_t ki = TMath::Nint(meani-3); + Int_t ki = TMath::Nint(meani); if (ki<0) ki=0; if (ki>=fMaxPad) ki = fMaxPad-1; - Int_t kj = TMath::Nint(meanj-3); + Int_t kj = TMath::Nint(meanj); if (kj<0) kj=0; if (kj>=fMaxTime-3) kj=fMaxTime-4; - // ki and kj shifted to "real" coordinata + // ki and kj shifted as integers coordinata if (fRowDig) { c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0); c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1); c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2); } - + c.SetRow(fRow); + c.SetDetector(fSector); Float_t s2 = c.GetSigmaY2(); Float_t w=fParam->GetPadPitchWidth(fSector); - c.SetSigmaY2(s2*w*w); s2 = c.GetSigmaZ2(); - w=fZWidth; - c.SetSigmaZ2(s2*w*w); - c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector)); - c.SetZ(fZWidth*(meanj-3)); - c.SetZ(c.GetZ() - 3.*fParam->GetZSigma()); // PASA delay - c.SetZ(fSign*(fParam->GetZLength() - c.GetZ())); + c.SetSigmaZ2(s2*fZWidth*fZWidth); + // + // + // + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + if (!transform) { + AliFatal("Tranformations not in calibDB"); + } + transform->SetCurrentRecoParam((AliTPCRecoParam*)fRecoParam); + Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()}; + Int_t i[1]={fSector}; + transform->Transform(x,i,0,1); + c.SetX(x[0]); + c.SetY(x[1]); + c.SetZ(x[2]); + // + // if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) { - //c.SetSigmaY2(c.GetSigmaY2()*25.); - //c.SetSigmaZ2(c.GetSigmaZ2()*4.); c.SetType(-(c.GetType()+3)); //edge clusters } if (fLoop==2) c.SetType(100); + if (!AcceptCluster(&c)) return; - TClonesArray * arr = fRowCl->GetArray(); - // AliTPCclusterMI * cl = - new ((*arr)[fNcluster]) AliTPCclusterMI(c); + // select output + TClonesArray * arr = 0; + AliTPCclusterMI * cl = 0; + + if(fBClonesArray==kFALSE) { + arr = fRowCl->GetArray(); + cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c); + } else { + cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c); + } + + // if (fRecoParam->DumpSignal() &&matrix ) { +// Int_t nbins=0; +// Float_t *graph =0; +// if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin() &&fBDumpSignal){ +// nbins = fMaxTime; +// graph = &(fBins[fMaxTime*(pos/fMaxTime)]); +// } +// AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph); +// cl->SetInfo(info); +// } + if (!fRecoParam->DumpSignal()) { + cl->SetInfo(0); + } + + if (AliTPCReconstructor::StreamLevel()>1) { + Float_t xyz[3]; + cl->GetGlobalXYZ(xyz); + (*fDebugStreamer)<<"Clusters"<< + "Cl.="<5) { + AliInfo("Parameter Dumps"); + fParam->Dump(); + fRecoParam->Dump(); } + AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); AliSimDigits digarr, *dummy=&digarr; fRowDig = dummy; fInput->GetBranch("Segment")->SetAddress(&dummy); Stat_t nentries = fInput->GetEntries(); - fMaxTime=fParam->GetMaxTBin()+6; // add 3 virtual time bins before and 3 after + fMaxTime=fRecoParam->GetLastBin()+6; // add 3 virtual time bins before and 3 after Int_t nclusters = 0; - + for (Int_t n=0; nGetEvent(n); - Int_t row; - if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,row)) { + if (!fParam->AdjustSectorRow(digarr.GetID(),fSector,fRow)) { cerr<<"AliTPC warning: invalid segment ID ! "<SetClass("AliTPCclusterMI"); - clrow->SetArray(1); + Int_t row = fRow; + AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector + AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector + // - clrow->SetID(digarr.GetID()); - fOutput->GetBranch("Segment")->SetAddress(&clrow); + fRowCl->SetID(digarr.GetID()); + if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl); fRx=fParam->GetPadRowRadii(fSector,row); @@ -504,241 +752,784 @@ void AliTPCclustererMI::Digits2Clusters() fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after - fBins =new Int_t[fMaxBin]; - fResBins =new Int_t[fMaxBin]; //fBins with residuals after 1 finder loop - memset(fBins,0,sizeof(Int_t)*fMaxBin); + fBins =new Float_t[fMaxBin]; + fSigBins =new Int_t[fMaxBin]; + fNSigBins = 0; + memset(fBins,0,sizeof(Float_t)*fMaxBin); if (digarr.First()) //MI change do { - Short_t dig=digarr.CurrentDigit(); + Float_t dig=digarr.CurrentDigit(); if (dig<=fParam->GetZeroSup()) continue; Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3; - fBins[i*fMaxTime+j]=dig; + Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn()); + Int_t bin = i*fMaxTime+j; + if (gain>0){ + fBins[bin]=dig/gain; + }else{ + fBins[bin]=0; + } + fSigBins[fNSigBins++]=bin; } while (digarr.Next()); digarr.ExpandTrackBuffer(); - FindClusters(); - - fOutput->Fill(); - delete clrow; + FindClusters(noiseROC); + FillRow(); + fRowCl->GetArray()->Clear("C"); nclusters+=fNcluster; - delete[] fBins; - delete[] fResBins; - } + delete[] fBins; + delete[] fSigBins; + } + Info("Digits2Clusters", "Number of found clusters : %d", nclusters); } +void AliTPCclustererMI::ProcessSectorData(){ + // + // Process the data for the current sector + // + + AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); + AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector + AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector + //check the presence of the calibration + if (!noiseROC ||!pedestalROC ) { + AliError(Form("Missing calibration per sector\t%d\n",fSector)); + return; + } + Int_t nRows=fParam->GetNRow(fSector); + Bool_t calcPedestal = fRecoParam->GetCalcPedestal(); + Int_t zeroSup = fParam->GetZeroSup(); + // if (calcPedestal) { + if (kFALSE ) { + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad = fParam->GetNPads(fSector, iRow); + + for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) { + // + // Temporary fix for data production - !!!! MARIAN + // The noise calibration should take mean and RMS - currently the Gaussian fit used + // In case of double peak - the pad should be rejected + // + // Line mean - if more than given digits over threshold - make a noise calculation + // and pedestal substration + if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue; + // + if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data + Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3]; + //Float_t pedestal = TMath::Median(fMaxTime, p); + Int_t id[3] = {fSector, iRow, iPad-3}; + // calib values + Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3); + Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3); + Double_t rmsEvent = rmsCalib; + Double_t pedestalEvent = pedestalCalib; + ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent); + if (rmsEventGetFirstBin()) + fAllBins[iRow][bin] = 0; + if (iTimeBin > fRecoParam->GetLastBin()) + fAllBins[iRow][bin] = 0; + if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup) + fAllBins[iRow][bin] = 0; + if (fAllBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS + fAllBins[iRow][bin] = 0; + if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin; + } + } + } + } + + if (AliTPCReconstructor::StreamLevel()>5) { + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad = fParam->GetNPads(fSector,iRow); + + for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) { + for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) { + Int_t bin = iPad*fMaxTime+iTimeBin; + Float_t signal = fAllBins[iRow][bin]; + if (AliTPCReconstructor::StreamLevel()>3 && signal>3) { + Double_t x[]={iRow,iPad-3,iTimeBin-3}; + Int_t i[]={fSector}; + AliTPCTransform trafo; + trafo.Transform(x,i,0,1); + Double_t gx[3]={x[0],x[1],x[2]}; + trafo.RotatedGlobal2Global(fSector,gx); + // fAllSigBins[iRow][fAllNSigBins[iRow]++] + Int_t rowsigBins = fAllNSigBins[iRow]; + Int_t first=fAllSigBins[iRow][0]; + Int_t last= 0; + // if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1]; + + if (AliTPCReconstructor::StreamLevel()>5) { + (*fDebugStreamer)<<"Digits"<< + "sec="<SetID(fParam->GetIndex(fSector, fRow)); + if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl); + + fRx = fParam->GetPadRowRadii(fSector, fRow); + fPadLength = fParam->GetPadPitchLength(fSector, fRow); + fPadWidth = fParam->GetPadPitchWidth(); + fMaxPad = fParam->GetNPads(fSector,fRow); + fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after + + fBins = fAllBins[fRow]; + fSigBins = fAllSigBins[fRow]; + fNSigBins = fAllNSigBins[fRow]; + + FindClusters(noiseROC); + + FillRow(); + if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C"); + fNclusters += fNcluster; + + } // End of loop to find clusters +} + + void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) { //----------------------------------------------------------------- -// This is a cluster finder for raw data. +// This is a cluster finder for the TPC raw data. +// The method assumes NO ordering of the altro channels. +// The pedestal subtraction can be switched on and off +// using an option of the TPC reconstructor //----------------------------------------------------------------- + fRecoParam = AliTPCReconstructor::GetRecoParam(); + if (!fRecoParam){ + AliFatal("Can not get the reconstruction parameters"); + } + if(AliTPCReconstructor::StreamLevel()>5) { + AliInfo("Parameter Dumps"); + fParam->Dump(); + fRecoParam->Dump(); + } + fRowDig = NULL; - if (!fOutput) { - Error("Digits2Clusters", "output tree not initialised"); - return; + AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); + AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); + // + AliTPCRawStreamV3 input(rawReader,(AliAltroMapping**)mapping); + fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader(); + if (fEventHeader){ + fTimeStamp = fEventHeader->Get("Timestamp"); + fEventType = fEventHeader->Get("Type"); + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + transform->SetCurrentTimeStamp(fTimeStamp); + transform->SetCurrentRun(rawReader->GetRunNumber()); + } + + // creaate one TClonesArray for all clusters + if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000); + // reset counter + fNclusters = 0; + + fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after +// const Int_t kNIS = fParam->GetNInnerSector(); +// const Int_t kNOS = fParam->GetNOuterSector(); +// const Int_t kNS = kNIS + kNOS; + fZWidth = fParam->GetZWidth(); + Int_t zeroSup = fParam->GetZeroSup(); + // + // Clean-up + // + AliTPCROC * roc = AliTPCROC::Instance(); + Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1); + Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1); + for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { + // + Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + fAllNSigBins[iRow]=0; } + Int_t prevSector=-1; rawReader->Reset(); - AliTPCRawStream input(rawReader); + Int_t digCounter=0; + // + // Loop over DDLs + // + const Int_t kNIS = fParam->GetNInnerSector(); + const Int_t kNOS = fParam->GetNOuterSector(); + const Int_t kNS = kNIS + kNOS; + + for(fSector = 0; fSector < kNS; fSector++) { + + Int_t nRows = 0; + Int_t nDDLs = 0, indexDDL = 0; + if (fSector < kNIS) { + nRows = fParam->GetNRowLow(); + fSign = (fSector < kNIS/2) ? 1 : -1; + nDDLs = 2; + indexDDL = fSector * 2; + } + else { + nRows = fParam->GetNRowUp(); + fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; + nDDLs = 4; + indexDDL = (fSector-kNIS) * 4 + kNIS * 2; + } + + // load the raw data for corresponding DDLs + rawReader->Reset(); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + + while (input.NextDDL()){ + if (input.GetSector() != fSector) + AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector())); + + //Int_t nRows = fParam->GetNRow(fSector); + + AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector + // Begin loop over altro data + Bool_t calcPedestal = fRecoParam->GetCalcPedestal(); + Float_t gain =1; + + //loop over pads + while ( input.NextChannel() ) { + Int_t iRow = input.GetRow(); + if (iRow < 0){ + continue; + } + if (iRow >= nRows){ + AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !", + iRow, 0, nRows -1)); + continue; + } + //pad + Int_t iPad = input.GetPad(); + if (iPad < 0 || iPad >= nPadsMax) { + AliError(Form("Pad index (%d) outside the range (%d -> %d) !", + iPad, 0, nPadsMax-1)); + continue; + } + gain = gainROC->GetValue(iRow,iPad); + iPad+=3; + + //loop over bunches + while ( input.NextBunch() ){ + Int_t startTbin = (Int_t)input.GetStartTimeBin(); + Int_t bunchlength = (Int_t)input.GetBunchLength(); + const UShort_t *sig = input.GetSignals(); + for (Int_t iTime = 0; iTimeGetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){ + continue; + AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !", + iTimeBin, 0, iTimeBin -1)); + } + iTimeBin+=3; + //signal + Float_t signal=(Float_t)sig[iTime]; + if (!calcPedestal && signal <= zeroSup) continue; + + if (!calcPedestal) { + Int_t bin = iPad*fMaxTime+iTimeBin; + if (gain>0){ + fAllBins[iRow][bin] = signal/gain; + }else{ + fAllBins[iRow][bin] =0; + } + fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin; + }else{ + fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal; + } + fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal + + // Temporary + digCounter++; + }// end loop signals in bunch + }// end loop bunches + } // end loop pads + // + // + // + // + // Now loop over rows and perform pedestal subtraction + if (digCounter==0) continue; + } // End of loop over sectors + //process last sector + if ( digCounter>0 ){ + ProcessSectorData(); + for (Int_t iRow = 0; iRow < fParam->GetNRow(fSector); iRow++) { + Int_t maxPad = fParam->GetNPads(fSector,iRow); + Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + fAllNSigBins[iRow] = 0; + } + prevSector=fSector; + digCounter=0; + } + } + + if (rawReader->GetEventId() && fOutput ){ + Info("Digits2Clusters", "File %s Event\t%u\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters); + } + + if(rawReader->GetEventId()) { + Info("Digits2Clusters", "Event\t%u\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters); + } + + if(fBClonesArray) { + //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast()); + } +} + + + + +void AliTPCclustererMI::Digits2ClustersOld +(AliRawReader* rawReader) +{ +//----------------------------------------------------------------- +// This is a cluster finder for the TPC raw data. +// The method assumes NO ordering of the altro channels. +// The pedestal subtraction can be switched on and off +// using an option of the TPC reconstructor +//----------------------------------------------------------------- + fRecoParam = AliTPCReconstructor::GetRecoParam(); + if (!fRecoParam){ + AliFatal("Can not get the reconstruction parameters"); + } + if(AliTPCReconstructor::StreamLevel()>5) { + AliInfo("Parameter Dumps"); + fParam->Dump(); + fRecoParam->Dump(); + } fRowDig = NULL; - Int_t nclusters = 0; + AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); + AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); + // + AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping); + fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader(); + if (fEventHeader){ + fTimeStamp = fEventHeader->Get("Timestamp"); + fEventType = fEventHeader->Get("Type"); + } + + // creaate one TClonesArray for all clusters + if(fBClonesArray && !fOutputClonesArray) fOutputClonesArray = new TClonesArray("AliTPCclusterMI",1000); + // reset counter + fNclusters = 0; - fMaxTime = fParam->GetMaxTBin() + 6; // add 3 virtual time bins before and 3 after + fMaxTime = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after const Int_t kNIS = fParam->GetNInnerSector(); const Int_t kNOS = fParam->GetNOuterSector(); const Int_t kNS = kNIS + kNOS; fZWidth = fParam->GetZWidth(); Int_t zeroSup = fParam->GetZeroSup(); + // + // Clean-up + // - fBins = NULL; - Int_t** splitRows = new Int_t* [kNS*2]; - Int_t** splitRowsRes = new Int_t* [kNS*2]; - for (Int_t iSector = 0; iSector < kNS*2; iSector++) - splitRows[iSector] = NULL; - Int_t iSplitRow = -1; - - Bool_t next = kTRUE; - while (next) { - next = input.Next(); - - // when the sector or row number has changed ... - if (input.IsNewRow() || !next) { - - // ... find clusters in the previous pad row, and ... - if (fBins) { - if ((iSplitRow < 0) || splitRows[fSector + kNS*iSplitRow]) { - fRowCl = new AliTPCClustersRow; - fRowCl->SetClass("AliTPCclusterMI"); - fRowCl->SetArray(1); - fRowCl->SetID(fParam->GetIndex(fSector, input.GetPrevRow())); - fOutput->GetBranch("Segment")->SetAddress(&fRowCl); - - FindClusters(); - - fOutput->Fill(); - delete fRowCl; - nclusters += fNcluster; - delete[] fBins; - delete[] fResBins; - if (iSplitRow >= 0) splitRows[fSector + kNS*iSplitRow] = NULL; - - } else if (iSplitRow >= 0) { - splitRows[fSector + kNS*iSplitRow] = fBins; - splitRowsRes[fSector + kNS*iSplitRow] = fResBins; - } - } + AliTPCROC * roc = AliTPCROC::Instance(); + Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1); + Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1); + for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { + // + Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + fAllNSigBins[iRow]=0; + } + // + // Loop over sectors + // + for(fSector = 0; fSector < kNS; fSector++) { - if (!next) break; + Int_t nRows = 0; + Int_t nDDLs = 0, indexDDL = 0; + if (fSector < kNIS) { + nRows = fParam->GetNRowLow(); + fSign = (fSector < kNIS/2) ? 1 : -1; + nDDLs = 2; + indexDDL = fSector * 2; + } + else { + nRows = fParam->GetNRowUp(); + fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; + nDDLs = 4; + indexDDL = (fSector-kNIS) * 4 + kNIS * 2; + } - // ... prepare for the next pad row - fSector = input.GetSector(); - Int_t iRow = input.GetRow(); - fRx = fParam->GetPadRowRadii(fSector, iRow); + // load the raw data for corresponding DDLs + rawReader->Reset(); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + + // select only good sector + input.Next(); + if(input.GetSector() != fSector) continue; + + AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector - iSplitRow = -1; - if (fSector < kNIS) { - fMaxPad = fParam->GetNPadsLow(iRow); - fSign = (fSector < kNIS/2) ? 1 : -1; - if (iRow == 30) iSplitRow = 0; - } else { - fMaxPad = fParam->GetNPadsUp(iRow); - fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; - if (iRow == 27) iSplitRow = 0; - else if (iRow == 76) iSplitRow = 1; - } - fPadLength = fParam->GetPadPitchLength(fSector, iRow); - fPadWidth = fParam->GetPadPitchWidth(); + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad; + if (fSector < kNIS) + maxPad = fParam->GetNPadsLow(iRow); + else + maxPad = fParam->GetNPadsUp(iRow); + + Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + fAllNSigBins[iRow] = 0; + } - fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after - if ((iSplitRow < 0) || !splitRows[fSector + kNS*iSplitRow]) { - fBins = new Int_t[fMaxBin]; - fResBins = new Int_t[fMaxBin]; //fBins with residuals after 1 finder loop - memset(fBins, 0, sizeof(Int_t)*fMaxBin); - } else { - fBins = splitRows[fSector + kNS*iSplitRow]; - fResBins = splitRowsRes[fSector + kNS*iSplitRow]; + Int_t digCounter=0; + // Begin loop over altro data + Bool_t calcPedestal = fRecoParam->GetCalcPedestal(); + Float_t gain =1; + Int_t lastPad=-1; + + input.Reset(); + while (input.Next()) { + if (input.GetSector() != fSector) + AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector())); + + + Int_t iRow = input.GetRow(); + if (iRow < 0){ + continue; } - } - // fill fBins with digits data - if (input.GetSignal() <= zeroSup) continue; - Int_t i = input.GetPad() + 3; - Int_t j = input.GetTime() + 3; - fBins[i*fMaxTime+j] = input.GetSignal(); - } - - // find clusters in split rows that were skipped until now. - // this can happen if the rows were not splitted - for (fSector = 0; fSector < kNS; fSector++) - for (Int_t iSplit = 0; iSplit < 2; iSplit++) - if (splitRows[fSector + kNS*iSplit]) { - - Int_t iRow = -1; - if (fSector < kNIS) { - iRow = 30; - fMaxPad = fParam->GetNPadsLow(iRow); - fSign = (fSector < kNIS/2) ? 1 : -1; - } else { - if (iSplit == 0) iRow = 27; else iRow = 76; - fMaxPad = fParam->GetNPadsUp(iRow); - fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; - } - fRx = fParam->GetPadRowRadii(fSector, iRow); - fPadLength = fParam->GetPadPitchLength(fSector, iRow); - fPadWidth = fParam->GetPadPitchWidth(); - - fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after - fBins = splitRows[fSector + kNS*iSplit]; - fResBins = splitRowsRes[fSector + kNS*iSplit]; - - fRowCl = new AliTPCClustersRow; - fRowCl->SetClass("AliTPCclusterMI"); - fRowCl->SetArray(1); - fRowCl->SetID(fParam->GetIndex(fSector, iRow)); - fOutput->GetBranch("Segment")->SetAddress(&fRowCl); - - FindClusters(); - - fOutput->Fill(); - delete fRowCl; - nclusters += fNcluster; - delete[] fBins; - delete[] fResBins; + if (iRow < 0 || iRow >= nRows){ + AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !", + iRow, 0, nRows -1)); + continue; + } + //pad + Int_t iPad = input.GetPad(); + if (iPad < 0 || iPad >= nPadsMax) { + AliError(Form("Pad index (%d) outside the range (%d -> %d) !", + iPad, 0, nPadsMax-1)); + continue; + } + if (iPad!=lastPad){ + gain = gainROC->GetValue(iRow,iPad); + lastPad = iPad; + } + iPad+=3; + //time + Int_t iTimeBin = input.GetTime(); + if ( iTimeBin < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){ + continue; + AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !", + iTimeBin, 0, iTimeBin -1)); + } + iTimeBin+=3; + + //signal + Float_t signal = input.GetSignal(); + if (!calcPedestal && signal <= zeroSup) continue; + + if (!calcPedestal) { + Int_t bin = iPad*fMaxTime+iTimeBin; + if (gain>0){ + fAllBins[iRow][bin] = signal/gain; + }else{ + fAllBins[iRow][bin] =0; + } + fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin; + }else{ + fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal; } + fAllBins[iRow][iPad*fMaxTime+0]+=1.; // pad with signal + + // Temporary + digCounter++; + } // End of the loop over altro data + // + // + // + // + // Now loop over rows and perform pedestal subtraction + if (digCounter==0) continue; + ProcessSectorData(); + } // End of loop over sectors + + if (rawReader->GetEventId() && fOutput ){ + Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), fNclusters); + } + + if(rawReader->GetEventId()) { + Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), fNclusters); + } - delete[] splitRows; - delete[] splitRowsRes; - Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters); + if(fBClonesArray) { + //Info("Digits2Clusters", "Number of found clusters : %d\n",fOutputClonesArray->GetEntriesFast()); + } } -void AliTPCclustererMI::FindClusters() +void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC) { - //add virtual charge at the edge - for (Int_t i=0; i0){ - Float_t amp2 = fBins[i+4*fMaxTime]; - if (amp2==0) amp2=0.5; - Float_t sigma2 = GetSigmaY2(i); - amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); - if (gDebug>4) printf("\n%f\n",amp0); - } - fBins[i+2*fMaxTime] = Int_t(amp0); - amp0 = 0; - amp1 = fBins[(fMaxPad+2)*fMaxTime+i]; - if (amp1>0){ - Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime]; - if (amp2==0) amp2=0.5; - Float_t sigma2 = GetSigmaY2(i); - amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); - if (gDebug>4) printf("\n%f\n",amp0); - } - fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0); - } - -// memcpy(fResBins,fBins, fMaxBin*2); - memcpy(fResBins,fBins, fMaxBin); + // + // add virtual charge at the edge + // + Double_t kMaxDumpSize = 500000; + if (!fOutput) { + fBDumpSignal =kFALSE; + }else{ + if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag + } + fNcluster=0; - //first loop - for "gold cluster" fLoop=1; - Int_t *b=&fBins[-1]+2*fMaxTime; - Int_t crtime = Int_t((fParam->GetZLength()-AliTPCReconstructor::GetCtgRange()*fRx)/fZWidth-5); - - for (Int_t i=2*fMaxTime; iGetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5); + Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs(); + Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs(); + Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs(); + Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma(); + Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma(); + Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma(); + Int_t useOnePadCluster = fRecoParam->GetUseOnePadCluster(); + for (Int_t iSig = 0; iSig < fNSigBins; iSig++) { + Int_t i = fSigBins[iSig]; + if (i%fMaxTime<=crtime) continue; + Float_t *b = &fBins[i]; + //absolute custs + if (b[0]GetValue(fRow, i/fMaxTime); + if (noise>fRecoParam->GetMaxNoise()) continue; + // sigma cuts + if (b[0]GetMax()<400) return kTRUE; + Double_t ratio = cl->GetQ()/cl->GetMax(); + if (cl->GetMax()>700){ + if ((ratio - int(ratio)>0.8)) return kFALSE; } - */ + if ((ratio - int(ratio)<0.95)) return kTRUE; + return kFALSE; } + + +Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t id[3], Double_t &rmsEvent, Double_t &pedestalEvent){ + // + // process signal on given pad - + streaming of additional information in special mode + // + // id[0] - sector + // id[1] - row + // id[2] - pad + + // + // ESTIMATE pedestal and the noise + // + const Int_t kPedMax = 100; + Float_t max = 0; + Float_t maxPos = 0; + Int_t median = -1; + Int_t count0 = 0; + Int_t count1 = 0; + Float_t rmsCalib = rmsEvent; // backup initial value ( from calib) + Float_t pedestalCalib = pedestalEvent;// backup initial value ( from calib) + Int_t firstBin = fRecoParam->GetFirstBin(); + // + UShort_t histo[kPedMax]; + //memset(histo,0,kPedMax*sizeof(UShort_t)); + for (Int_t i=0; imax && i>firstBin) { + max = signal[i]; + maxPos = i; + } + if (signal[i]>kPedMax-1) continue; + histo[int(signal[i]+0.5)]++; + count0++; + } + // + for (Int_t i=1; ikPedMax) continue; + if (count06<0.6*count1){ + count06+=histo[median-idelta]; + mean06 +=histo[median-idelta]*(median-idelta); + rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta); + count06+=histo[median+idelta]; + mean06 +=histo[median+idelta]*(median+idelta); + rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta); + } + if (count09<0.9*count1){ + count09+=histo[median-idelta]; + mean09 +=histo[median-idelta]*(median-idelta); + rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta); + count09+=histo[median+idelta]; + mean09 +=histo[median+idelta]*(median+idelta); + rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta); + } + if (count10<0.95*count1){ + count10+=histo[median-idelta]; + mean +=histo[median-idelta]*(median-idelta); + rms +=histo[median-idelta]*(median-idelta)*(median-idelta); + count10+=histo[median+idelta]; + mean +=histo[median+idelta]*(median+idelta); + rms +=histo[median+idelta]*(median+idelta)*(median+idelta); + } + } + if (count10) { + mean /=count10; + rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean)); + } + if (count06) { + mean06/=count06; + rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06)); + } + if (count09) { + mean09/=count09; + rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09)); + } + rmsEvent = rms09; + // + pedestalEvent = median; + if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median; + // + UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])}; + // + // Dump mean signal info + // + if (AliTPCReconstructor::StreamLevel()>0) { + (*fDebugStreamer)<<"Signal"<< + "TimeStamp="<GetDumpAmplitudeMin(); // minimal signal to be dumped + Float_t *dsignal = new Float_t[nchannels]; + Float_t *dtime = new Float_t[nchannels]; + for (Int_t i=0; i0) { + if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin()) + (*fDebugStreamer)<<"SignalB"<< // pads with signal + "TimeStamp="<fRecoParam->GetMaxNoise()) { + pedestalEvent+=1024.; + return 1024+median; // sign noisy channel in debug mode + } + return median; +} + + + + +