X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCclustererMI.cxx;h=ecfad37cd71f0df786b9d237dfbf9459afa8680f;hb=49519b04d13fb4961e50bf1c6e63b8614e54f635;hp=e6e54bb8c933420b57ca8667e6dbb7a524d8c6fa;hpb=7fef31a60454688a0d0330cfbd12ebe5c68a44d4;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCclustererMI.cxx b/TPC/AliTPCclustererMI.cxx index e6e54bb8c93..ecfad37cd71 100644 --- a/TPC/AliTPCclustererMI.cxx +++ b/TPC/AliTPCclustererMI.cxx @@ -18,47 +18,117 @@ //------------------------------------------------------- // 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); +// 1.c HLT clusters - Digits2Clusters and Digits2Clusters(AliRawReader* rawReader) +// invoke ReadHLTClusters() +// +// fUseHLTClusters - switches between different inputs +// 1 -> only TPC raw/sim data +// 2 -> if present TPC raw/sim data, otherwise HLT clusters +// 3 -> only HLT clusters +// 4 -> if present HLT clusters, otherwise TPC raw/sim data +// +// 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 "TGraph.h" -#include "TF1.h" -#include "TRandom.h" -#include "AliMathBase.h" +#include +#include +#include +#include +#include +#include +#include +#include "TSystem.h" +#include "TClass.h" -#include "AliTPCClustersArray.h" -#include "AliTPCClustersRow.h" #include "AliDigits.h" -#include "AliSimDigits.h" -#include "AliTPCParam.h" -#include "AliTPCRecoParam.h" +#include "AliLoader.h" +#include "AliLog.h" +#include "AliMathBase.h" +#include "AliRawEventHeaderBase.h" #include "AliRawReader.h" -#include "AliTPCRawStream.h" #include "AliRunLoader.h" -#include "AliLoader.h" -#include "Riostream.h" -#include -#include "AliTPCcalibDB.h" +#include "AliSimDigits.h" #include "AliTPCCalPad.h" #include "AliTPCCalROC.h" -#include "TTreeStream.h" -#include "AliLog.h" - +#include "AliTPCClustersArray.h" +#include "AliTPCClustersRow.h" +#include "AliTPCParam.h" +#include "AliTPCRawStream.h" +#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, const AliTPCRecoParam * recoParam) +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), + fUseHLTClusters(1), + fAllBins(NULL), + fAllSigBins(NULL), + fAllNSigBins(NULL) { - fIsOldRCUFormat = kFALSE; + // + // COSNTRUCTOR + // param - tpc parameters for given file + // recoparam - reconstruction parameters + // fInput =0; - fOutput=0; fParam = par; if (recoParam) { fRecoParam = recoParam; @@ -67,14 +137,116 @@ ClassImp(AliTPCclustererMI) fRecoParam = AliTPCReconstructor::GetRecoParam(); if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam(); } - fDebugStreamer = new TTreeSRedirector("TPCsignal.root"); - fAmplitudeHisto = 0; -} + + 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), + fUseHLTClusters(1), + 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(){ - DumpHistos(); - if (fAmplitudeHisto) delete fAmplitudeHisto; + // + // + // 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) @@ -93,16 +265,32 @@ 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; + if (!tree) return; + fOutput= tree; + AliTPCClustersRow clrow("AliTPCclusterMI"); AliTPCClustersRow *pclrow=&clrow; - clrow.SetClass("AliTPCclusterMI"); - clrow.SetArray(1); // to make Clones array 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()+fParam->GetNTBinsL1()*fParam->GetZWidth(); @@ -118,7 +306,7 @@ Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){ //sigma z2 = in digits - angle estimated supposing vertex constraint 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; @@ -130,16 +318,20 @@ Float_t AliTPCclustererMI::GetSigmaZ2(Int_t iz){ 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}; Float_t * matrix[5]; //5x5 matrix with digits - indexing i = 0 ..4 j = -2..2 - Float_t * resmatrix[5]; 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); @@ -216,7 +408,7 @@ AliTPCclusterMI &c) // if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return; - if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) { + if ( ((ry <1.2) && (rz<1.2)) || (!fRecoParam->GetDoUnfold())) { // //if cluster looks like expected or Unfolding not switched on //standard COG is used @@ -227,18 +419,12 @@ 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] -= 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; } // @@ -262,18 +448,16 @@ 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]); } @@ -310,7 +494,7 @@ void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[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]; @@ -401,8 +585,6 @@ void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5] } } } - if (gDebug>4) - printf("%f\n", recmatrix[2][2]); } @@ -425,53 +607,95 @@ 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() + fParam->GetNTBinsL1()*fParam->GetZWidth()); // PASA delay + L1 delay - c.SetZ(fSign*(fParam->GetZLength() - c.GetZ())); - c.SetX(fRx); - c.SetDetector(fSector); - c.SetRow(fRow); - + c.SetSigmaZ2(s2*fZWidth*fZWidth); + // + // + // + + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + if (!transform) { + AliFatal("Tranformations not in calibDB"); + return; + } + 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; + + // 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); + } - TClonesArray * arr = fRowCl->GetArray(); - // 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); + } + + 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(); + //----------------------------------------------------------------- + // Use HLT clusters + //----------------------------------------------------------------- + fUseHLTClusters = fRecoParam->GetUseHLTClusters(); + AliInfo(Form("Usage of HLT clusters in TPC reconstruction : %d",fUseHLTClusters)); + + if (fUseHLTClusters == 3 || fUseHLTClusters == 4) { + AliInfo("Using HLT clusters for TPC off-line reconstruction"); + fZWidth = fParam->GetZWidth(); + Int_t iResult = ReadHLTClusters(); + + // HLT clusters present + if (iResult >= 0 && fNclusters > 0) + return; + + // HLT clusters not present + if (iResult < 0 || fNclusters == 0) { + if (fUseHLTClusters == 3) { + AliError("No HLT clusters present, but requested."); + return; + } + else { + AliInfo("Now trying to read TPC RAW"); + } + } + // Some other problem during cluster reading + else { + AliWarning("Some problem while unpacking of HLT clusters."); + return; + } + } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) { + + //----------------------------------------------------------------- + // Run TPC off-line clusterer + //----------------------------------------------------------------- + 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; @@ -513,15 +777,11 @@ void AliTPCclustererMI::Digits2Clusters() } Int_t row = fRow; AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector - + AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector // - AliTPCClustersRow *clrow= new AliTPCClustersRow(); - fRowCl = clrow; - clrow->SetClass("AliTPCclusterMI"); - clrow->SetArray(1); - 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); @@ -529,7 +789,7 @@ void AliTPCclustererMI::Digits2Clusters() fZWidth = fParam->GetZWidth(); if (fSector < kNIS) { fMaxPad = fParam->GetNPadsLow(row); - fSign = (fSector < kNIS/2) ? 1 : -1; + fSign = (fSector < kNIS/2) ? 1 : -1; fPadLength = fParam->GetPadPitchLength(fSector,row); fPadWidth = fParam->GetPadPitchWidth(); } else { @@ -542,9 +802,9 @@ void AliTPCclustererMI::Digits2Clusters() fMaxBin=fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after fBins =new Float_t[fMaxBin]; - fResBins =new Float_t[fMaxBin]; //fBins with residuals after 1 finder loop + fSigBins =new Int_t[fMaxBin]; + fNSigBins = 0; memset(fBins,0,sizeof(Float_t)*fMaxBin); - memset(fResBins,0,sizeof(Float_t)*fMaxBin); if (digarr.First()) //MI change do { @@ -552,22 +812,169 @@ void AliTPCclustererMI::Digits2Clusters() if (dig<=fParam->GetZeroSup()) continue; Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3; Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn()); - fBins[i*fMaxTime+j]=dig/gain; + 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); + + if (fUseHLTClusters == 2 && nclusters == 0) { + AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters."); + + fZWidth = fParam->GetZWidth(); + ReadHLTClusters(); + } } +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) { //----------------------------------------------------------------- @@ -576,33 +983,292 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) // 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; + //----------------------------------------------------------------- + // Use HLT clusters + //----------------------------------------------------------------- + fUseHLTClusters = fRecoParam->GetUseHLTClusters(); + + AliInfo(Form("Usage of HLT clusters in TPC reconstruction : %d",fUseHLTClusters)); + + if (fUseHLTClusters == 3 || fUseHLTClusters == 4) { + AliInfo("Using HLT clusters for TPC off-line reconstruction"); + fZWidth = fParam->GetZWidth(); + Int_t iResult = ReadHLTClusters(); + + // HLT clusters present + if (iResult >= 0 && fNclusters > 0) + return; + + // HLT clusters not present + if (iResult < 0 || fNclusters == 0) { + if (fUseHLTClusters == 3) { + AliError("No HLT clusters present, but requested."); + return; + } + else { + AliInfo("Now trying to read TPC RAW"); + } + } + // Some other problem during cluster reading + else { + AliWarning("Some problem while unpacking of HLT clusters."); + return; + } + } // if (fUseHLTClusters == 3 || fUseHLTClusters == 4) { + + //----------------------------------------------------------------- + // Run TPC off-line clusterer + //----------------------------------------------------------------- + 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(); + 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()); + } + + if (fUseHLTClusters == 2 && fNclusters == 0) { + AliInfo("No clusters from TPC Raw data, now trying to read HLT clusters."); + + fZWidth = fParam->GetZWidth(); + ReadHLTClusters(); + } +} + + + + +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; 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"); + } - Int_t nclusters = 0; + // 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 + // - Float_t** allBins = NULL; - Float_t** allBinsRes = NULL; - + 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++) { - AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector - Int_t nRows = 0; Int_t nDDLs = 0, indexDDL = 0; if (fSector < kNIS) { @@ -618,321 +1284,283 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) indexDDL = (fSector-kNIS) * 4 + kNIS * 2; } - allBins = new Float_t*[nRows]; - allBinsRes = new Float_t*[nRows]; + // load the raw data for corresponding DDLs + rawReader->Reset(); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + // select only good sector + if (!input.Next()) continue; + if(input.GetSector() != fSector) continue; + + AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector + for (Int_t iRow = 0; iRow < nRows; iRow++) { Int_t maxPad; if (fSector < kNIS) - maxPad = fParam->GetNPadsLow(iRow); + maxPad = fParam->GetNPadsLow(iRow); else - maxPad = fParam->GetNPadsUp(iRow); + maxPad = fParam->GetNPadsUp(iRow); Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after - allBins[iRow] = new Float_t[maxBin]; - allBinsRes[iRow] = new Float_t[maxBin]; - memset(allBins[iRow],0,sizeof(Float_t)*maxBin); - memset(allBinsRes[iRow],0,sizeof(Float_t)*maxBin); + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + fAllNSigBins[iRow] = 0; } - // Loas the raw data for corresponding DDLs - rawReader->Reset(); - AliTPCRawStream input(rawReader); - input.SetOldRCUFormat(fIsOldRCUFormat); - rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); - + Int_t digCounter=0; // Begin loop over altro data - while (input.Next()) { + 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 || iRow >= nRows) - AliFatal(Form("Pad-row index (%d) outside the range (%d -> %d) !", - iRow, 0, nRows -1)); - - Int_t iPad = input.GetPad() + 3; - - Int_t maxPad; - if (fSector < kNIS) - maxPad = fParam->GetNPadsLow(iRow); - else - maxPad = fParam->GetNPadsUp(iRow); - - if (input.GetPad() < 0 || input.GetPad() >= maxPad) - AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !", - input.GetPad(), 0, maxPad -1)); + AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector())); - Int_t iTimeBin = input.GetTime() + 3; - if ( input.GetTime() < 0 || input.GetTime() >= fParam->GetMaxTBin()) - AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !", - input.GetTime(), 0, fParam->GetMaxTBin() -1)); - Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after + Int_t iRow = input.GetRow(); + if (iRow < 0){ + continue; + } - if (((iPad*fMaxTime+iTimeBin) >= maxBin) || - ((iPad*fMaxTime+iTimeBin) < 0)) - AliFatal(Form("Index outside the allowed range" - " Sector=%d Row=%d Pad=%d Timebin=%d" - " (Max.index=%d)",fSector,iRow,iPad,iTimeBin,maxBin)); + 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 (!fPedSubtraction && signal <= zeroSup) continue; - if (!fRecoParam->GetCalcPedestal() && signal <= zeroSup) continue; + 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 - Float_t gain = gainROC->GetValue(iRow,input.GetPad()); - allBins[iRow][iPad*fMaxTime+iTimeBin] = signal/gain; - 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 (fPedSubtraction) { - if (fRecoParam->GetCalcPedestal()) { - 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 = 0; iPad < maxPad + 6; iPad++) { - 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}; - Float_t pedestal = ProcesSignal(p, fMaxTime, id); - for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) { - allBins[iRow][iPad*fMaxTime+iTimeBin] -= pedestal; - if (iTimeBin < AliTPCReconstructor::GetRecoParam()->GetFirstBin()) - allBins[iRow][iPad*fMaxTime+iTimeBin] = 0; - if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin()) - allBins[iRow][iPad*fMaxTime+iTimeBin] = 0; - if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup) - allBins[iRow][iPad*fMaxTime+iTimeBin] = 0; - } - } - } - } - - // 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)); - 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]; - fResBins = allBinsRes[fRow]; - - FindClusters(); - - fOutput->Fill(); - delete fRowCl; - nclusters += fNcluster; - } // End of loop to find clusters - - - for (Int_t iRow = 0; iRow < nRows; iRow++) { - delete [] allBins[iRow]; - delete [] allBinsRes[iRow]; - } - - delete [] allBins; - delete [] allBinsRes; - + if (digCounter==0) continue; + ProcessSectorData(); } // End of loop over sectors - Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters); + 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); + } + 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] = 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] = 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; - Float_t *b=&fBins[-1]+2*fMaxTime; - Int_t crtime = Int_t((fParam->GetZLength()-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-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 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 - Int_t offset =100; - Float_t kMin =fRecoParam->GetDumpAmplitudeMin(); // minimal signal to be dumped - Double_t median = TMath::Median(nchannels-offset, &(signal[offset])); - if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median; - // + // id[2] - pad - Double_t mean = TMath::Mean(nchannels-offset, &(signal[offset])); - Double_t rms = TMath::RMS(nchannels-offset, &(signal[offset])); - Double_t *dsignal = new Double_t[nchannels]; - Double_t *dtime = new Double_t[nchannels]; - Float_t max = 0; - Float_t maxPos = 0; - for (Int_t i=0; imax && i 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++; } - - Double_t mean06=0, mean09=0; - Double_t rms06=0, rms09=0; - AliMathBase::EvaluateUni(nchannels-offset, &(dsignal[offset]), mean06, rms06, int(0.6*(nchannels-offset))); - AliMathBase::EvaluateUni(nchannels-offset, &(dsignal[offset]), mean09, rms09, int(0.9*(nchannels-offset))); // + for (Int_t i=1; iGetNSectors() - && uid[1]< AliTPCROC::Instance()->GetNRows(uid[0]) && - uid[2] < AliTPCROC::Instance()->GetNPads(uid[0], uid[1])){ - if (!fAmplitudeHisto){ - fAmplitudeHisto = new TObjArray(72); - } - TObjArray * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]); - if (!sectorArray){ - Int_t npads = AliTPCROC::Instance()->GetNChannels(uid[0]); - sectorArray = new TObjArray(npads); - fAmplitudeHisto->AddAt(sectorArray, uid[0]); + Float_t count10=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ; + Float_t count06=histo[median] ,mean06=histo[median]*median, rms06=histo[median]*median*median ; + Float_t count09=histo[median] ,mean09=histo[median]*median, rms09=histo[median]*median*median ; + // + for (Int_t idelta=1; idelta<10; idelta++){ + if (median-idelta<=0) continue; + if (median+idelta>kPedMax) 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); } - Int_t position = uid[2]+ AliTPCROC::Instance()->GetRowIndexes(uid[0])[uid[1]]; - TH1F * histo = (TH1F*)sectorArray->UncheckedAt(position); - if (!histo){ - char hname[100]; - sprintf(hname,"Amp_%d_%d_%d",uid[0],uid[1],uid[2]); - TFile * backup = gFile; - fDebugStreamer->GetFile()->cd(); - histo = new TH1F(hname, hname, 100, 5,100); - //histo->SetDirectory(0); // histogram not connected to directory -(File) - sectorArray->AddAt(histo, position); - if (backup) backup->cd(); + 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); } - for (Int_t i=0; i0) histo->Fill(signal[i]); + 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); } } - // - TGraph * graph; - Bool_t random = (gRandom->Rndm()<0.0001); - if (max-median>kMin || rms06>2.*fParam->GetZeroSup() || random){ - graph =new TGraph(nchannels, dtime, dsignal); - if (rms06>2.*fParam->GetZeroSup() || random) - (*fDebugStreamer)<<"SignalN"<< //noise pads - or random sample of pads - "Sector="<kMin) - (*fDebugStreamer)<<"SignalB"<< // pads with signal - "Sector="<0) { + (*fDebugStreamer)<<"Signal"<< + "TimeStamp="<ceThreshold) for (Int_t i=nchannels-2; i>1; i--){ - if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){ - cemaxpos=i; - break; - } - } - if (cemaxpos!=0){ - for (Int_t i=cemaxpos-cemin; i0 && iceQmax) ceQmax=val; - } - } - if (ceQmax&&ceQsum>ceSumThreshold) { - ceTime/=ceQsum; - (*fDebugStreamer)<<"Signalce"<< - "Sector="<ggThreshold) for (Int_t i=1; iggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){ - ggmaxpos=i; - break; - } + Float_t kMin =fRecoParam->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 && i0){ - Double_t val=dsignal[i]- ggmean; - ggTime+=val*dtime[i]; - ggQsum+=val; - if (val>ggQmax) ggQmax=val; - } - } - if (ggQmax&&ggQsum>ggSumThreshold) { - ggTime/=ggQsum; - (*fDebugStreamer)<<"Signalgg"<< - "Sector="<0) { + if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin()) + (*fDebugStreamer)<<"SignalB"<< // pads with signal + "TimeStamp="<fRecoParam->GetMaxNoise()) return 1024+median; // sign noisy channel in debug mode + if (rms06>fRecoParam->GetMaxNoise()) { + pedestalEvent+=1024.; + return 1024+median; // sign noisy channel in debug mode + } return median; } +Int_t AliTPCclustererMI::ReadHLTClusters() +{ + // + // read HLT clusters instead of off line custers, + // used in Digits2Clusters + // + TObject* pClusterAccess=NULL; + TClass* pCl=NULL; + ROOT::NewFunc_t pNewFunc=NULL; + do { + pCl=TClass::GetClass("AliHLTTPCClusterAccessHLTOUT"); + } while (!pCl && gSystem->Load("libAliHLTTPC.so")==0); + if (!pCl || (pNewFunc=pCl->GetNew())==NULL) { + AliError("can not load class description of AliHLTTPCClusterAccessHLTOUT, aborting ..."); + return -1; + } + + void* p=(*pNewFunc)(NULL); + if (!p) { + AliError("unable to create instance of AliHLTTPCClusterAccessHLTOUT"); + return -2; + } + pClusterAccess=reinterpret_cast(p); + if (!pClusterAccess) { + AliError("instance not of type TObject"); + return -3 ; + } -void AliTPCclustererMI::DumpHistos(){ - if (!fAmplitudeHisto) return; - for (UInt_t isector=0; isectorGetNSectors(); isector++){ - TObjArray * array = (TObjArray*)fAmplitudeHisto->UncheckedAt(isector); - if (!array) continue; - for (UInt_t ipad = 0; ipad <(UInt_t)array->GetEntriesFast(); ipad++){ - TH1F * histo = (TH1F*) array->UncheckedAt(ipad); - if (!histo) continue; - if (histo->GetEntries()<100) continue; - histo->Fit("gaus","q"); - Float_t mean = histo->GetMean(); - Float_t rms = histo->GetRMS(); - Float_t gmean = histo->GetFunction("gaus")->GetParameter(1); - Float_t gsigma = histo->GetFunction("gaus")->GetParameter(2); - Float_t max = histo->GetFunction("gaus")->GetParameter(0); - - // get pad number - UInt_t row=0, pad =0; - const UInt_t *indexes = AliTPCROC::Instance()->GetRowIndexes(isector); - for (UInt_t irow=0; irow< AliTPCROC::Instance()->GetNRows(isector); irow++){ - if (indexes[irow]GetNPads(isector,row))/2; - // - (*fDebugStreamer)<<"Fit"<< - "Sector="<UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad)); + const Int_t kNIS = fParam->GetNInnerSector(); + const Int_t kNOS = fParam->GetNOuterSector(); + const Int_t kNS = kNIS + kNOS; + fNclusters = 0; + + for(fSector = 0; fSector < kNS; fSector++) { + + Int_t iResult = 1; + TString param("sector="); param+=fSector; + pClusterAccess->Clear(); + pClusterAccess->Execute("read", param, &iResult); + if (iResult < 0) { + return iResult; + AliError("HLT Clusters can not be found"); } - } + + if (pClusterAccess->FindObject("clusterarray")==NULL) { + AliError("HLT clusters requested, but not cluster array not present"); + return -4; + } + + TClonesArray* clusterArray=dynamic_cast(pClusterAccess->FindObject("clusterarray")); + if (!clusterArray) { + AliError("HLT cluster array is not of class type TClonesArray"); + return -5; + } + + AliDebug(4,Form("Reading %d clusters from HLT for sector %d", clusterArray->GetEntriesFast(), fSector)); + + Int_t nClusterSector=0; + Int_t nRows=fParam->GetNRow(fSector); + + for (fRow = 0; fRow < nRows; fRow++) { + fRowCl->SetID(fParam->GetIndex(fSector, fRow)); + if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl); + fNcluster=0; // reset clusters per row + + 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]; + + for (Int_t i=0; iGetEntriesFast(); i++) { + if (!clusterArray->At(i)) + continue; + + AliTPCclusterMI* cluster=dynamic_cast(clusterArray->At(i)); + if (!cluster) continue; + if (cluster->GetRow()!=fRow) continue; + nClusterSector++; + AddCluster(*cluster, NULL, 0); + } + + FillRow(); + fRowCl->GetArray()->Clear("c"); + + } // for (fRow = 0; fRow < nRows; fRow++) { + if (nClusterSector!=clusterArray->GetEntriesFast()) { + AliError(Form("Failed to read %d out of %d HLT clusters", + clusterArray->GetEntriesFast()-nClusterSector, + clusterArray->GetEntriesFast())); + } + fNclusters+=nClusterSector; + } // for(fSector = 0; fSector < kNS; fSector++) { + + delete pClusterAccess; + + Info("Digits2Clusters", "Number of converted HLT clusters : %d", fNclusters); + + return 0; }