X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCclustererMI.cxx;h=ecfad37cd71f0df786b9d237dfbf9459afa8680f;hb=19dcf5042148839eaec09a5a00f17c17c11c615b;hp=ed26f3c405212663ee916affa9b72c13cfdec7bd;hpb=a1ec4d073607f1e7bad0e027746dad890100bcc6;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCclustererMI.cxx b/TPC/AliTPCclustererMI.cxx index ed26f3c4052..ecfad37cd71 100644 --- a/TPC/AliTPCclustererMI.cxx +++ b/TPC/AliTPCclustererMI.cxx @@ -18,6 +18,31 @@ //------------------------------------------------------- // 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 //------------------------------------------------------- @@ -27,10 +52,12 @@ #include #include #include +#include #include #include #include -#include +#include "TSystem.h" +#include "TClass.h" #include "AliDigits.h" #include "AliLoader.h" @@ -46,11 +73,13 @@ #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) @@ -59,10 +88,11 @@ ClassImp(AliTPCclustererMI) AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam): fBins(0), - fResBins(0), + fSigBins(0), + fNSigBins(0), fLoop(0), fMaxBin(0), - fMaxTime(0), + fMaxTime(1006), // 1000>940 so use 1000, add 3 virtual time bins before and 3 after fMaxPad(0), fSector(-1), fRow(-1), @@ -72,30 +102,33 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar fPadLength(0), fZWidth(0), fPedSubtraction(kFALSE), - fIsOldRCUFormat(kFALSE), fEventHeader(0), fTimeStamp(0), fEventType(0), fInput(0), fOutput(0), + fOutputArray(0), + fOutputClonesArray(0), fRowCl(0), fRowDig(0), fParam(0), fNcluster(0), - fAmplitudeHisto(0), + fNclusters(0), fDebugStreamer(0), fRecoParam(0), fBDumpSignal(kFALSE), - fFFTr2c(0) + fBClonesArray(kFALSE), + fUseHLTClusters(1), + fAllBins(NULL), + fAllSigBins(NULL), + fAllNSigBins(NULL) { // // COSNTRUCTOR // param - tpc parameters for given file // recoparam - reconstruction parameters // - fIsOldRCUFormat = kFALSE; fInput =0; - fOutput=0; fParam = par; if (recoParam) { fRecoParam = recoParam; @@ -104,16 +137,40 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar fRecoParam = AliTPCReconstructor::GetRecoParam(); if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam(); } - fDebugStreamer = new TTreeSRedirector("TPCsignal.root"); - fAmplitudeHisto = 0; - Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin(); - fFFTr2c = TVirtualFFT::FFT(1, &nPoints, "R2C K"); + + 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), - fResBins(0), + fSigBins(0), + fNSigBins(0), fLoop(0), fMaxBin(0), fMaxTime(0), @@ -126,19 +183,26 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI ¶m) fPadLength(0), fZWidth(0), fPedSubtraction(kFALSE), - fIsOldRCUFormat(kFALSE), fEventHeader(0), fTimeStamp(0), fEventType(0), fInput(0), fOutput(0), + fOutputArray(0), + fOutputClonesArray(0), fRowCl(0), fRowDig(0), fParam(0), fNcluster(0), - fAmplitudeHisto(0), + fNclusters(0), fDebugStreamer(0), - fRecoParam(0) + fRecoParam(0), + fBDumpSignal(kFALSE), + fBClonesArray(kFALSE), + fUseHLTClusters(1), + fAllBins(NULL), + fAllSigBins(NULL), + fAllNSigBins(NULL) { // // dummy @@ -156,9 +220,33 @@ AliTPCclustererMI & AliTPCclustererMI::operator =(const AliTPCclustererMI & para } //______________________________________________________________ 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) @@ -177,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(); @@ -226,10 +330,8 @@ AliTPCclusterMI &c) // 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); @@ -306,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 @@ -317,21 +419,12 @@ AliTPCclusterMI &c) meanj +=j0; //set cluster parameters c.SetQ(sumw); - c.SetY(meani*fPadWidth); - c.SetZ(meanj*fZWidth); - c.SetPad(meani); - c.SetTimeBin(meanj); + c.SetPad(meani-2.5); + c.SetTimeBin(meanj-3); c.SetSigmaY2(mi2); c.SetSigmaZ2(mj2); + c.SetType(0); AddCluster(c,(Float_t*)vmatrix,k); - //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; - return; } // @@ -355,10 +448,8 @@ AliTPCclusterMI &c) meanj +=j0; //set cluster parameters c.SetQ(sumu); - c.SetY(meani*fPadWidth); - c.SetZ(meanj*fZWidth); - c.SetPad(meani); - c.SetTimeBin(meanj); + c.SetPad(meani-2.5); + c.SetTimeBin(meanj-3); c.SetSigmaY2(mi2); c.SetSigmaZ2(mj2); c.SetType(Char_t(overlap)+1); @@ -367,8 +458,6 @@ AliTPCclusterMI &c) //unfolding 2 meani-=i0; meanj-=j0; - if (gDebug>4) - printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]); } @@ -405,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]; @@ -496,8 +585,6 @@ void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5] } } } - if (gDebug>4) - printf("%f\n", recmatrix[2][2]); } @@ -520,66 +607,94 @@ Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, F return max; } -void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t pos){ +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)); - if (!fRecoParam->GetBYMirror()){ - if (fSector%36>17){ - c.SetY(-(meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector)); - } + c.SetSigmaZ2(s2*fZWidth*fZWidth); + // + // + // + + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; + if (!transform) { + AliFatal("Tranformations not in calibDB"); + return; } - 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(fSector) - c.GetZ())); - c.SetX(fRx); - c.SetDetector(fSector); - c.SetRow(fRow); - + 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); - if (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); + // 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(); } + //----------------------------------------------------------------- + // 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; @@ -625,13 +779,9 @@ void AliTPCclustererMI::Digits2Clusters() 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); @@ -639,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 { @@ -652,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 { @@ -662,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(noiseROC); - - fOutput->Fill(); - delete clrow; + 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) { //----------------------------------------------------------------- @@ -686,58 +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; - AliTPCROC * roc = AliTPCROC::Instance(); + AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); - AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals(); - AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); - AliTPCRawStream input(rawReader); + 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(); // - //alocate memory for sector - maximal case + // 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); - allBins = new Float_t*[nRowsMax]; - allBinsRes = new Float_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]; - allBinsRes[iRow] = new Float_t[maxBin]; - memset(allBins[iRow],0,sizeof(Float_t)*maxBin); + 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 - AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector - AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector - Int_t nRows = 0; Int_t nDDLs = 0, indexDDL = 0; if (fSector < kNIS) { @@ -753,141 +1284,111 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) indexDDL = (fSector-kNIS) * 4 + kNIS * 2; } + // 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 - memset(allBins[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(); - 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; + + input.Reset(); while (input.Next()) { - digCounter++; if (input.GetSector() != fSector) - AliFatal(Form("Sector index mismatch ! Expected (%d), but got (%d) !",fSector,input.GetSector())); + 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) !", + if (iRow < 0){ + continue; + } + + 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) - AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !", + 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; + gain = gainROC->GetValue(iRow,iPad); + lastPad = iPad; } iPad+=3; //time Int_t iTimeBin = input.GetTime(); - if ( iTimeBin < 0 || iTimeBin >= fParam->GetMaxTBin()) - AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !", + 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 && signal <= zeroSup) continue; + if (!calcPedestal) { - allBins[iRow][iPad*fMaxTime+iTimeBin] = signal/gain; + 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{ - allBins[iRow][iPad*fMaxTime+iTimeBin] = signal; + fAllBins[iRow][iPad*fMaxTime+iTimeBin] = signal; } - allBins[iRow][iPad*fMaxTime+0]=1.; // pad with 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; - // if (fPedSubtraction) { - if (calcPedestal) { - 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++) { - 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][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; - if (allBins[iRow][iPad*fMaxTime+iTimeBin] < 3.0*rmsEvent) // 3 sigma cut on RMS - 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(noiseROC); - - fOutput->Fill(); - delete fRowCl; - nclusters += fNcluster; - } // End of loop to find clusters + 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); + } - for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { - delete [] allBins[iRow]; - delete [] allBinsRes[iRow]; - } - delete [] allBins; - delete [] allBinsRes; - - Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters); + 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(AliTPCCalROC * noiseROC) @@ -897,37 +1398,14 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC) // add virtual charge at the edge // Double_t kMaxDumpSize = 500000; - if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag - // - if (0) 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; + if (!fOutput) { + fBDumpSignal =kFALSE; + }else{ + if (fRecoParam->GetCalcPedestal() && fOutput->GetZipBytes()< kMaxDumpSize) fBDumpSignal =kTRUE; //dump signal flag } - memcpy(fResBins,fBins, fMaxBin*sizeof(Float_t)); - // - // - // + fNcluster=0; fLoop=1; - Float_t *b=&fBins[-1]+2*fMaxTime; Int_t crtime = Int_t((fParam->GetZLength(fSector)-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5); Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs(); Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs(); @@ -935,32 +1413,32 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC) Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma(); Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma(); Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma(); - for (Int_t i=2*fMaxTime; iGetUseOnePadCluster(); + 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){ // @@ -981,7 +1475,6 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t // ESTIMATE pedestal and the noise // const Int_t kPedMax = 100; - Double_t kMaxDebugSize = 5000000.; Float_t max = 0; Float_t maxPos = 0; Int_t median = -1; @@ -989,10 +1482,11 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t 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 = AliTPCReconstructor::GetRecoParam()->GetFirstBin(); + Int_t firstBin = fRecoParam->GetFirstBin(); // UShort_t histo[kPedMax]; - memset(histo,0,kPedMax*sizeof(UShort_t)); + //memset(histo,0,kPedMax*sizeof(UShort_t)); + for (Int_t i=0; imax && i>firstBin) { @@ -1042,12 +1536,18 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t rms +=histo[median+idelta]*(median+idelta)*(median+idelta); } } - mean /=count10; - mean06/=count06; - mean09/=count09; - rms = TMath::Sqrt(TMath::Abs(rms/count10-mean*mean)); - rms06 = TMath::Sqrt(TMath::Abs(rms06/count06-mean06*mean06)); - rms09 = TMath::Sqrt(TMath::Abs(rms09/count09-mean09*mean09)); + 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; @@ -1057,7 +1557,8 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t // // Dump mean signal info // - (*fDebugStreamer)<<"Signal"<< + if (AliTPCReconstructor::StreamLevel()>0) { + (*fDebugStreamer)<<"Signal"<< "TimeStamp="<GetNSectors() - && uid[1]< roc->GetNRows(uid[0]) && - uid[2] GetNPads(uid[0], uid[1])){ - TObjArray * sectorArray = (TObjArray*)fAmplitudeHisto->UncheckedAt(uid[0]); - if (!sectorArray){ - Int_t npads =roc->GetNChannels(uid[0]); - sectorArray = new TObjArray(npads); - fAmplitudeHisto->AddAt(sectorArray, uid[0]); - } - Int_t position = uid[2]+roc->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(); - } - for (Int_t i=0; iFill(signal[i]); - } - } // // // @@ -1119,265 +1591,31 @@ Double_t AliTPCclustererMI::ProcesSignal(Float_t *signal, Int_t nchannels, Int_t dtime[i] = i; dsignal[i] = signal[i]; } - // - // Digital noise - // - // if (max-median>30.*TMath::Max(1.,Double_t(rms06)) && (((*fDebugStreamer)<<"SignalDN").GetSize()GetFirstBin(); -// Int_t njumps1=0; -// Double_t deltaT1[2000]; -// Double_t deltaA1[2000]; -// Int_t lastJump1 = fRecoParam->GetFirstBin(); -// Int_t njumps2=0; -// Double_t deltaT2[2000]; -// Double_t deltaA2[2000]; -// Int_t lastJump2 = fRecoParam->GetFirstBin(); - -// for (Int_t itime=fRecoParam->GetFirstBin()+1; itimeGetLastBin()-1; itime++){ -// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) && -// TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) && -// (dsignal[itime-1]-median<5.*rms06) && -// (dsignal[itime+1]-median<5.*rms06) -// ){ -// deltaA0[njumps0] = dsignal[itime]-dsignal[itime-1]; -// deltaT0[njumps0] = itime-lastJump0; -// lastJump0 = itime; -// njumps0++; -// } -// if (TMath::Abs(dsignal[itime]-dsignal[itime-1])>30.*TMath::Max(1.,Double_t(rms06)) && -// (dsignal[itime-1]-median<5.*rms06) -// ) { -// deltaA1[njumps1] = dsignal[itime]-dsignal[itime-1]; -// deltaT1[njumps1] = itime-lastJump1; -// lastJump1 = itime; -// njumps1++; -// } -// if (TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,Double_t(rms06)) && -// (dsignal[itime+1]-median<5.*rms06) -// ) { -// deltaA2[njumps2] = dsignal[itime]-dsignal[itime+1]; -// deltaT2[njumps2] = itime-lastJump2; -// lastJump2 = itime; -// njumps2++; -// } -// } -// // -// if (njumps0>0 || njumps1>0 || njumps2>0){ -// TGraph *graphDN0 = new TGraph(njumps0, deltaT0, deltaA0); -// TGraph *graphDN1 = new TGraph(njumps1, deltaT1, deltaA1); -// TGraph *graphDN2 = new TGraph(njumps2, deltaT2, deltaA2); -// (*fDebugStreamer)<<"SignalDN"<< //digital - noise pads - or random sample of pads -// "TimeStamp="<Rndm()<0.0003); - if (((*fDebugStreamer)<<"SignalN").GetSize()kMin || rms06>1.*fParam->GetZeroSup() || random){ - graph =new TGraph(nchannels, dtime, dsignal); - if (rms06>1.*fParam->GetZeroSup() || random){ - //Double_t *input, Double_t threshold, Bool_t locMax, Double_t *freq, Double_t *re, Double_t *im, Double_t *mag, Double_t *phi); - Float_t * input = &(dsignal[fRecoParam->GetFirstBin()]); - Float_t freq[2000], re[2000], im[2000], mag[2000], phi[2000]; - Int_t npoints = TransformFFT(input, -1,kFALSE, freq, re, im, mag, phi); - TGraph *graphMag0 = new TGraph(npoints, freq, mag); - TGraph *graphPhi0 = new TGraph(npoints, freq, phi); - npoints = TransformFFT(input, 0.5,kTRUE, freq, re, im, mag, phi); - TGraph *graphMag1 = new TGraph(npoints, freq, mag); - TGraph *graphPhi1 = new TGraph(npoints, freq, phi); - - (*fDebugStreamer)<<"SignalN"<< //noise pads - or random sample of pads - "TimeStamp="<kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin()) - (*fDebugStreamer)<<"SignalB"<< // pads with signal - "TimeStamp="<nchannels/2; i--){ - if ( (dsignal[i]-mean06)>ceThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] ){ - cemaxpos=i; - break; - } - } - if (cemaxpos!=0){ - ceQmax = 0; - Int_t cemaxpos2=0; - for (Int_t i=cemaxpos-20; inchannels-1) continue; - Double_t val=dsignal[i]- cemean; - if (val>ceQmax){ - cemaxpos2=i; - ceQmax = val; - } - } - cemaxpos = cemaxpos2; - - for (Int_t i=cemaxpos-kCemin; i0 && i0){ - Double_t val=dsignal[i]- cemean; - ceTime+=val*dtime[i]; - ceQsum+=val; - if (val>ceQmax) ceQmax=val; - } - } - if (ceQmax&&ceQsum>ceSumThreshold) { - ceTime/=ceQsum; - (*fDebugStreamer)<<"Signalce"<< - "TimeStamp="<ggThreshold && dsignal[i]>=dsignal[i+1] && dsignal[i]>=dsignal[i-1] && - (dsignal[i]+dsignal[i+1]+dsignal[i-1]-3*mean06)>ggSumThreshold){ - ggmaxpos=i; - if (dsignal[i-1]>dsignal[i+1]) ggmaxpos=i-1; - break; - } - } - if (ggmaxpos!=0){ - for (Int_t i=ggmaxpos-1; 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"<< - "TimeStamp="<0) { + if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin()) + (*fDebugStreamer)<<"SignalB"<< // pads with signal + "TimeStamp="<GetNSectors(); 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 gmeanErr = histo->GetFunction("gaus")->GetParError(1); - Float_t gsigmaErr = histo->GetFunction("gaus")->GetParError(2); - Float_t max = histo->GetFunction("gaus")->GetParameter(0); - - // get pad number - UInt_t row=0, pad =0; - const UInt_t *indexes =roc->GetRowIndexes(isector); - for (UInt_t irow=0; irowGetNRows(isector); irow++){ - if (indexes[irow]<=ipad){ - row = irow; - pad = ipad-indexes[irow]; - } - } - Int_t rpad = pad - (AliTPCROC::Instance()->GetNPads(isector,row))/2; - // - (*fDebugStreamer)<<"Fit"<< - "TimeStamp="<UncheckedAt(ipad)) fDebugStreamer->StoreObject(array->UncheckedAt(ipad)); - } + + 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 ; } -} + 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"); + } -Int_t AliTPCclustererMI::TransformFFT(Float_t *input, Float_t threshold, Bool_t locMax, Float_t *freq, Float_t *re, Float_t *im, Float_t *mag, Float_t *phi) -{ - // - // calculate fourrie transform - // return only frequncies with mag over threshold - // if locMax is spectified only freque with local maxima over theshold is returned - - if (! fFFTr2c) return kFALSE; - if (!freq) return kFALSE; - - Int_t current=0; - Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin(); - Double_t *in = new Double_t[nPoints]; - Double_t *rfft = new Double_t[nPoints]; - Double_t *ifft = new Double_t[nPoints]; - for (Int_t i=0; iSetPoints(in); - fFFTr2c->Transform(); - fFFTr2c->GetPointsComplex(rfft, ifft); - for (Int_t i=3; ilmag) continue; - if ( TMath::Sqrt(rfft[i+1]*rfft[i+1]+ifft[i+1]*ifft[i+1])/nPoints>lmag) continue; - if ( TMath::Sqrt(rfft[i-2]*rfft[i-2]+ifft[i-2]*ifft[i-2])/nPoints>lmag) continue; - if ( TMath::Sqrt(rfft[i+2]*rfft[i+2]+ifft[i+2]*ifft[i+2])/nPoints>lmag) continue; - if ( TMath::Sqrt(rfft[i-3]*rfft[i-3]+ifft[i-3]*ifft[i-3])/nPoints>lmag) continue; - if ( TMath::Sqrt(rfft[i+3]*rfft[i+3]+ifft[i+3]*ifft[i+3])/nPoints>lmag) continue; + if (pClusterAccess->FindObject("clusterarray")==NULL) { + AliError("HLT clusters requested, but not cluster array not present"); + return -4; } - - freq[current] = Float_t(i)/Float_t(nPoints); - // - re[current] = rfft[i]; - im[current] = ifft[i]; - mag[current]=lmag; - phi[current]=TMath::ATan2(ifft[i],rfft[i]); - current++; - } - delete [] in; - delete [] rfft; - delete [] ifft; - return current; -} + 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; +}