X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;ds=sidebyside;f=TPC%2FAliTPCclustererMI.cxx;h=6da558b7aa5f072a5c7c1feedbbe5bfa734fbe66;hb=e4042305eece1ea7ad68e8dda2b76a957e06df1e;hp=2e7cf06d148550c985c0af286377a1fe9d08605c;hpb=6c024a0e3ec23e58ab5ebb92466dec3bf801fbc9;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCclustererMI.cxx b/TPC/AliTPCclustererMI.cxx index 2e7cf06d148..6da558b7aa5 100644 --- a/TPC/AliTPCclustererMI.cxx +++ b/TPC/AliTPCclustererMI.cxx @@ -18,37 +18,169 @@ //------------------------------------------------------- // 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 +// +// 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 "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 "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(0), + fMaxPad(0), + fSector(-1), + fRow(-1), + fSign(0), + fRx(0), + fPadWidth(0), + fPadLength(0), + fZWidth(0), + fPedSubtraction(kFALSE), + fIsOldRCUFormat(kFALSE), + fEventHeader(0), + fTimeStamp(0), + fEventType(0), + fInput(0), + fOutput(0), + fOutputArray(0), + fRowCl(0), + fRowDig(0), + fParam(0), + fNcluster(0), + fDebugStreamer(0), + fRecoParam(0), + fBDumpSignal(kFALSE) { + // + // COSNTRUCTOR + // param - tpc parameters for given file + // recoparam - reconstruction parameters + // + fIsOldRCUFormat = kFALSE; 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(); + } + fDebugStreamer = new TTreeSRedirector("TPCsignal.root"); + Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin(); +} +//______________________________________________________________ +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), + fIsOldRCUFormat(kFALSE), + fEventHeader(0), + fTimeStamp(0), + fEventType(0), + fInput(0), + fOutput(0), + fOutputArray(0), + fRowCl(0), + fRowDig(0), + fParam(0), + fNcluster(0), + fDebugStreamer(0), + fRecoParam(0), + fBDumpSignal(kFALSE) +{ + // + // 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; + } +} + void AliTPCclustererMI::SetInput(TTree * tree) { // @@ -65,8 +197,11 @@ 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; + if (!tree) return; + fOutput= tree; AliTPCClustersRow clrow; AliTPCClustersRow *pclrow=&clrow; clrow.SetClass("AliTPCclusterMI"); @@ -75,9 +210,24 @@ void AliTPCclustererMI::SetOutput(TTree * tree) } +void AliTPCclustererMI::FillRow(){ + // + // fill the output container - + // 2 Options possible + // Tree + // TObjArray + // + if (fOutput) fOutput->Fill(); + if (!fOutput){ + // + if (!fOutputArray) fOutputArray = new TObjArray; + if (fRowCl) 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 +238,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 +249,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); @@ -188,8 +342,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 +353,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,12 +382,12 @@ 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; @@ -248,7 +398,7 @@ AliTPCclusterMI &c) -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 ) { // @@ -395,50 +545,77 @@ 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"); + } + 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 (!fRecoParam->GetBYMirror()){ + if (fSector%36>17){ + c.SetY(-c.GetY()); + } + } + 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); TClonesArray * arr = fRowCl->GetArray(); - // AliTPCclusterMI * cl = - new ((*arr)[fNcluster]) AliTPCclusterMI(c); + AliTPCclusterMI * cl = new ((*arr)[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); + } fNcluster++; } @@ -456,35 +633,33 @@ void AliTPCclustererMI::Digits2Clusters() return; } - if (!fOutput) { - Error("Digits2Clusters", "output tree not initialised"); - return; - } - + 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); - - clrow->SetID(digarr.GetID()); - fOutput->GetBranch("Segment")->SetAddress(&clrow); + Int_t row = fRow; + AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector + AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector + // + fRowCl= new AliTPCClustersRow(); + fRowCl->SetClass("AliTPCclusterMI"); + fRowCl->SetArray(1); + + fRowCl->SetID(digarr.GetID()); + if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl); fRx=fParam->GetPadRowRadii(fSector,row); @@ -504,26 +679,29 @@ 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; + fBins[bin]=dig/gain; + fSigBins[fNSigBins++]=bin; } while (digarr.Next()); digarr.ExpandTrackBuffer(); - FindClusters(); - - fOutput->Fill(); - delete clrow; + FindClusters(noiseROC); + FillRow(); + delete fRowCl; nclusters+=fNcluster; - delete[] fBins; - delete[] fResBins; + delete[] fBins; + delete[] fSigBins; } Info("Digits2Clusters", "Number of found clusters : %d", nclusters); @@ -532,213 +710,468 @@ void AliTPCclustererMI::Digits2Clusters() 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 //----------------------------------------------------------------- - if (!fOutput) { - Error("Digits2Clusters", "output tree not initialised"); - return; - } - - rawReader->Reset(); - AliTPCRawStream input(rawReader); fRowDig = NULL; + AliTPCROC * roc = AliTPCROC::Instance(); + AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); + AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); + AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); + // + AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping); + fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader(); + if (fEventHeader){ + fTimeStamp = fEventHeader->Get("Timestamp"); + fEventType = fEventHeader->Get("Type"); + } + Int_t nclusters = 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(); + // + //alocate memory for sector - maximal case + // + Float_t** allBins = NULL; + Int_t** allSigBins = NULL; + Int_t* allNSigBins = NULL; + Int_t nRowsMax = roc->GetNRows(roc->GetNSector()-1); + Int_t nPadsMax = roc->GetNPads(roc->GetNSector()-1,nRowsMax-1); + allBins = new Float_t*[nRowsMax]; + allSigBins = new Int_t*[nRowsMax]; + allNSigBins = 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 + allBins[iRow] = new Float_t[maxBin]; + memset(allBins[iRow],0,sizeof(Float_t)*maxBin); + allSigBins[iRow] = new Int_t[maxBin]; + allNSigBins[iRow]=0; + } + // + // Loop over sectors + // + for(fSector = 0; fSector < kNS; fSector++) { + + AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector + 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)); + continue; + } + 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; + } - 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; - } - } - - if (!next) break; - - // ... prepare for the next pad row - fSector = input.GetSector(); - Int_t iRow = input.GetRow(); - fRx = fParam->GetPadRowRadii(fSector, iRow); + 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(allBins[iRow],0,sizeof(Float_t)*maxBin); + allNSigBins[iRow] = 0; + } - 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; + // Loas the raw data for corresponding DDLs + rawReader->Reset(); + input.SetOldRCUFormat(fIsOldRCUFormat); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + Int_t digCounter=0; + // Begin loop over altro data + Bool_t calcPedestal = fRecoParam->GetCalcPedestal(); + Float_t gain =1; + Int_t lastPad=-1; + 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 || iRow >= nRows){ + AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !", + iRow, 0, nRows -1)); + continue; } - fPadLength = fParam->GetPadPitchLength(fSector, iRow); - fPadWidth = fParam->GetPadPitchWidth(); - - 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]; + //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; } - } - - // 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; + 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; + allBins[iRow][bin] = signal/gain; + allSigBins[iRow][allNSigBins[iRow]++] = bin; + }else{ + allBins[iRow][iPad*fMaxTime+iTimeBin] = signal; + } + allBins[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; + // if (calcPedestal) { + if (kTRUE) { + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad; + if (fSector < kNIS) + maxPad = fParam->GetNPadsLow(iRow); + else + maxPad = fParam->GetNPadsUp(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 && allBins[iRow][iPad*fMaxTime+0]<50) continue; + // + if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data + Float_t *p = &allBins[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()) + allBins[iRow][bin] = 0; + if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin()) + allBins[iRow][bin] = 0; + if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup) + allBins[iRow][bin] = 0; + if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS + allBins[iRow][bin] = 0; + if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin; + } } - 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; } + } + // Now loop over rows and find clusters + for (fRow = 0; fRow < nRows; fRow++) { + fRowCl = new AliTPCClustersRow; + fRowCl->SetClass("AliTPCclusterMI"); + fRowCl->SetArray(1); + fRowCl->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(); + if (fSector < kNIS) + fMaxPad = fParam->GetNPadsLow(fRow); + else + fMaxPad = fParam->GetNPadsUp(fRow); + fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after + + fBins = allBins[fRow]; + fSigBins = allSigBins[fRow]; + fNSigBins = allNSigBins[fRow]; - delete[] splitRows; - delete[] splitRowsRes; - Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters); + FindClusters(noiseROC); + FillRow(); + delete fRowCl; + nclusters += fNcluster; + } // End of loop to find clusters + } // End of loop over sectors + + for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { + delete [] allBins[iRow]; + delete [] allSigBins[iRow]; + } + delete [] allBins; + delete [] allSigBins; + delete [] allNSigBins; + +// if (rawReader->GetEventId() && fOutput ){ +// Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters); +// }else{ +// Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), nclusters); + +// } + } -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(); + 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]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 + // + (*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; ikMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin()) + (*fDebugStreamer)<<"SignalB"<< // pads with signal + "TimeStamp="<fRecoParam->GetMaxNoise()) { + pedestalEvent+=1024.; + return 1024+median; // sign noisy channel in debug mode + } + return median; } + + + + +