X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCclustererMI.cxx;h=5fa8076c484708ba758bc50b43ed8d5b64b0199a;hb=9e97e52de12f5ee493f484dd6aabc816071797b8;hp=a3cf16f101bcb21f39054e76389567ee23e3dcc9;hpb=9d4f75a9170a2cd09aa78a903c3d6a567485d5ea;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCclustererMI.cxx b/TPC/AliTPCclustererMI.cxx index a3cf16f101b..5fa8076c484 100644 --- a/TPC/AliTPCclustererMI.cxx +++ b/TPC/AliTPCclustererMI.cxx @@ -24,36 +24,141 @@ #include "AliTPCReconstructor.h" #include "AliTPCclustererMI.h" #include "AliTPCclusterMI.h" +#include "AliTPCclusterInfo.h" #include #include +#include "TGraph.h" +#include "TF1.h" +#include "TRandom.h" +#include "AliMathBase.h" + #include "AliTPCClustersArray.h" #include "AliTPCClustersRow.h" -#include "AliTPCRawStream.h" #include "AliDigits.h" #include "AliSimDigits.h" #include "AliTPCParam.h" +#include "AliTPCRecoParam.h" #include "AliRawReader.h" #include "AliTPCRawStream.h" +#include "AliRawEventHeaderBase.h" #include "AliRunLoader.h" #include "AliLoader.h" #include "Riostream.h" #include - #include "AliTPCcalibDB.h" #include "AliTPCCalPad.h" #include "AliTPCCalROC.h" - +#include "TTreeStream.h" +#include "AliLog.h" +#include "TVirtualFFT.h" ClassImp(AliTPCclustererMI) -AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par) +AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoParam * recoParam): + fBins(0), + fResBins(0), + fLoop(0), + fMaxBin(0), + fMaxTime(0), + fMaxPad(0), + fSector(-1), + fRow(-1), + fSign(0), + fRx(0), + fPadWidth(0), + fPadLength(0), + fZWidth(0), + fPedSubtraction(kFALSE), + fIsOldRCUFormat(kFALSE), + fEventHeader(0), + fTimeStamp(0), + fEventType(0), + fInput(0), + fOutput(0), + fRowCl(0), + fRowDig(0), + fParam(0), + fNcluster(0), + fAmplitudeHisto(0), + fDebugStreamer(0), + fRecoParam(0), + fFFTr2c(0) { + // + // COSNTRUCTOR + // param - tpc parameters for given file + // recoparam - reconstruction parameters + // + fIsOldRCUFormat = kFALSE; fInput =0; fOutput=0; fParam = par; + if (recoParam) { + fRecoParam = recoParam; + }else{ + //set default parameters if not specified + fRecoParam = AliTPCReconstructor::GetRecoParam(); + if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam(); + } + fDebugStreamer = new TTreeSRedirector("TPCsignal.root"); + fAmplitudeHisto = 0; + Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin(); + fFFTr2c = TVirtualFFT::FFT(1, &nPoints, "R2C K"); +} +//______________________________________________________________ +AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI ¶m) + :TObject(param), + fBins(0), + fResBins(0), + fLoop(0), + fMaxBin(0), + fMaxTime(0), + fMaxPad(0), + fSector(-1), + fRow(-1), + fSign(0), + fRx(0), + fPadWidth(0), + fPadLength(0), + fZWidth(0), + fPedSubtraction(kFALSE), + fIsOldRCUFormat(kFALSE), + fEventHeader(0), + fTimeStamp(0), + fEventType(0), + fInput(0), + fOutput(0), + fRowCl(0), + fRowDig(0), + fParam(0), + fNcluster(0), + fAmplitudeHisto(0), + fDebugStreamer(0), + fRecoParam(0) +{ + // + // 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; +} + void AliTPCclustererMI::SetInput(TTree * tree) { // @@ -107,6 +212,12 @@ 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 @@ -193,8 +304,10 @@ AliTPCclusterMI &c) // if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return; - if ( (ry <1.2) && (rz<1.2) ) { - //if cluster looks like expected + if ( (ry <1.2) && (rz<1.2) || (!fRecoParam->GetDoUnfold())) { + // + //if cluster looks like expected or Unfolding not switched on + //standard COG is used //+1.2 deviation from expected sigma accepted // c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2)); @@ -204,9 +317,11 @@ AliTPCclusterMI &c) c.SetQ(sumw); c.SetY(meani*fPadWidth); c.SetZ(meanj*fZWidth); + c.SetPad(meani); + c.SetTimeBin(meanj); c.SetSigmaY2(mi2); c.SetSigmaZ2(mj2); - AddCluster(c); + 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++){ @@ -214,6 +329,7 @@ AliTPCclusterMI &c) if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0; } resmatrix[2][0] =0; + return; } // @@ -239,10 +355,12 @@ AliTPCclusterMI &c) c.SetQ(sumu); c.SetY(meani*fPadWidth); c.SetZ(meanj*fZWidth); + c.SetPad(meani); + c.SetTimeBin(meanj); c.SetSigmaY2(mi2); c.SetSigmaZ2(mj2); c.SetType(Char_t(overlap)+1); - AddCluster(c); + AddCluster(c,(Float_t*)vmatrix,k); //unfolding 2 meani-=i0; @@ -400,7 +518,7 @@ 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 @@ -430,6 +548,11 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){ 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.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())); @@ -445,8 +568,17 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){ if (fLoop==2) c.SetType(100); TClonesArray * arr = fRowCl->GetArray(); - // AliTPCclusterMI * cl = - new ((*arr)[fNcluster]) AliTPCclusterMI(c); + AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c); + if (matrix) { + Int_t nbins=0; + Float_t *graph =0; + if (fRecoParam->GetCalcPedestal() && cl->GetMax()>fRecoParam->GetDumpAmplitudeMin()){ + nbins = fMaxTime; + graph = &(fBins[fMaxTime*(pos/fMaxTime)]); + } + AliTPCclusterInfo * info = new AliTPCclusterInfo(matrix,nbins,graph); + cl->SetInfo(info); + } fNcluster++; } @@ -470,6 +602,7 @@ void AliTPCclustererMI::Digits2Clusters() } AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); AliSimDigits digarr, *dummy=&digarr; fRowDig = dummy; @@ -488,7 +621,7 @@ 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; @@ -531,7 +664,7 @@ void AliTPCclustererMI::Digits2Clusters() } while (digarr.Next()); digarr.ExpandTrackBuffer(); - FindClusters(); + FindClusters(noiseROC); fOutput->Fill(); delete clrow; @@ -546,7 +679,10 @@ void AliTPCclustererMI::Digits2Clusters() void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) { //----------------------------------------------------------------- -// This is a cluster finder for raw data. +// This is a cluster finder for the TPC raw data. +// The method assumes NO ordering of the altro channels. +// The pedestal subtraction can be switched on and off +// using an option of the TPC reconstructor //----------------------------------------------------------------- if (!fOutput) { @@ -554,12 +690,18 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) return; } + fRowDig = NULL; + AliTPCROC * roc = AliTPCROC::Instance(); AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); - - rawReader->Reset(); + AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); AliTPCRawStream input(rawReader); + fEventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader(); + if (fEventHeader){ + fTimeStamp = fEventHeader->Get("Timestamp"); + fEventType = fEventHeader->Get("Type"); + } - fRowDig = NULL; Int_t nclusters = 0; @@ -569,136 +711,190 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) const Int_t kNS = kNIS + kNOS; fZWidth = fParam->GetZWidth(); Int_t zeroSup = fParam->GetZeroSup(); + // + //alocate memory for sector - maximal case + // + Float_t** allBins = NULL; + Float_t** allBinsRes = NULL; + 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); + } + // + // Loop over sectors + // + for(fSector = 0; fSector < kNS; fSector++) { - fBins = NULL; - Float_t** splitRows = new Float_t* [kNS*2]; - Float_t** splitRowsRes = new Float_t* [kNS*2]; - for (Int_t iSector = 0; iSector < kNS*2; iSector++) - splitRows[iSector] = NULL; - Int_t iSplitRow = -1; - - Bool_t next = kTRUE; - while (next) { - next = input.Next(); - - // when the sector or row number has changed ... - if (input.IsNewRow() || !next) { - - // ... find clusters in the previous pad row, and ... - if (fBins) { - if ((iSplitRow < 0) || splitRows[fSector + kNS*iSplitRow]) { - fRowCl = new AliTPCClustersRow; - fRowCl->SetClass("AliTPCclusterMI"); - fRowCl->SetArray(1); - fRowCl->SetID(fParam->GetIndex(fSector, input.GetPrevRow())); - fOutput->GetBranch("Segment")->SetAddress(&fRowCl); - - FindClusters(); - - fOutput->Fill(); - delete fRowCl; - nclusters += fNcluster; - delete[] fBins; - delete[] fResBins; - if (iSplitRow >= 0) splitRows[fSector + kNS*iSplitRow] = NULL; - - } else if (iSplitRow >= 0) { - splitRows[fSector + kNS*iSplitRow] = fBins; - splitRowsRes[fSector + kNS*iSplitRow] = fResBins; - } - } - - if (!next) break; + 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) { + nRows = fParam->GetNRowLow(); + fSign = (fSector < kNIS/2) ? 1 : -1; + nDDLs = 2; + indexDDL = fSector * 2; + } + else { + nRows = fParam->GetNRowUp(); + fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; + nDDLs = 4; + indexDDL = (fSector-kNIS) * 4 + kNIS * 2; + } - // ... prepare for the next pad row - fSector = input.GetSector(); - fRow = input.GetRow(); - Int_t iRow = input.GetRow(); - fRx = fParam->GetPadRowRadii(fSector, iRow); + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad; + if (fSector < kNIS) + maxPad = fParam->GetNPadsLow(iRow); + else + maxPad = fParam->GetNPadsUp(iRow); + + Int_t maxBin = fMaxTime*(maxPad+6); // add 3 virtual pads before and 3 after + memset(allBins[iRow],0,sizeof(Float_t)*maxBin); + } - iSplitRow = -1; - if (fSector < kNIS) { - fMaxPad = fParam->GetNPadsLow(iRow); - fSign = (fSector < kNIS/2) ? 1 : -1; - if (iRow == 30) iSplitRow = 0; - } else { - fMaxPad = fParam->GetNPadsUp(iRow); - fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; - if (iRow == 27) iSplitRow = 0; - else if (iRow == 76) iSplitRow = 1; + // Loas the raw data for corresponding DDLs + rawReader->Reset(); + input.SetOldRCUFormat(fIsOldRCUFormat); + rawReader->Select("TPC",indexDDL,indexDDL+nDDLs-1); + Int_t digCounter=0; + // Begin loop over altro data + Bool_t calcPedestal = fRecoParam->GetCalcPedestal(); + Float_t gain =1; + Int_t lastPad=-1; + while (input.Next()) { + digCounter++; + 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)); + //pad + Int_t iPad = input.GetPad(); + if (iPad < 0 || iPad >= nPadsMax) + AliFatal(Form("Pad index (%d) outside the range (%d -> %d) !", + iPad, 0, nPadsMax-1)); + if (iPad!=lastPad){ + gain = gainROC->GetValue(iRow,iPad); + lastPad = iPad; } - fPadLength = fParam->GetPadPitchLength(fSector, iRow); - fPadWidth = fParam->GetPadPitchWidth(); - - fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after - if ((iSplitRow < 0) || !splitRows[fSector + kNS*iSplitRow]) { - fBins = new Float_t[fMaxBin]; - fResBins = new Float_t[fMaxBin]; //fBins with residuals after 1 finder loop - memset(fBins, 0, sizeof(Float_t)*fMaxBin); - memset(fResBins, 0, sizeof(Float_t)*fMaxBin); - } else { - fBins = splitRows[fSector + kNS*iSplitRow]; - fResBins = splitRowsRes[fSector + kNS*iSplitRow]; + iPad+=3; + //time + Int_t iTimeBin = input.GetTime(); + if ( iTimeBin < 0 || iTimeBin >= fParam->GetMaxTBin()) + 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) { + allBins[iRow][iPad*fMaxTime+iTimeBin] = signal/gain; + }else{ + allBins[iRow][iPad*fMaxTime+iTimeBin] = signal; } - } - - // fill fBins with digits data - if (input.GetSignal() <= zeroSup) continue; - Int_t i = input.GetPad() + 3; - Int_t j = input.GetTime() + 3; - AliTPCCalROC * gainROC = gainTPC->GetCalROC(fSector); // pad gains per given sector - Float_t gain = gainROC->GetValue(fRow,input.GetPad()); - fBins[i*fMaxTime+j] = input.GetSignal()/gain; - } - - // find clusters in split rows that were skipped until now. - // this can happen if the rows were not splitted - for (fSector = 0; fSector < kNS; fSector++) - for (Int_t iSplit = 0; iSplit < 2; iSplit++) - if (splitRows[fSector + kNS*iSplit]) { - - Int_t iRow = -1; - if (fSector < kNIS) { - iRow = 30; - fMaxPad = fParam->GetNPadsLow(iRow); - fSign = (fSector < kNIS/2) ? 1 : -1; - } else { - if (iSplit == 0) iRow = 27; else iRow = 76; - fMaxPad = fParam->GetNPadsUp(iRow); - fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1; + allBins[iRow][iPad*fMaxTime+0]=1.; // pad with signal + } // 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; + } } - fRx = fParam->GetPadRowRadii(fSector, iRow); - fPadLength = fParam->GetPadPitchLength(fSector, iRow); - fPadWidth = fParam->GetPadPitchWidth(); - - fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after - fBins = splitRows[fSector + kNS*iSplit]; - fResBins = splitRowsRes[fSector + kNS*iSplit]; - - fRowCl = new AliTPCClustersRow; - fRowCl->SetClass("AliTPCclusterMI"); - fRowCl->SetArray(1); - fRowCl->SetID(fParam->GetIndex(fSector, iRow)); - fOutput->GetBranch("Segment")->SetAddress(&fRowCl); - - FindClusters(); - - fOutput->Fill(); - delete fRowCl; - nclusters += fNcluster; - delete[] fBins; - delete[] fResBins; } + } + // Now loop over rows and find clusters + for (fRow = 0; fRow < nRows; fRow++) { + fRowCl = new AliTPCClustersRow; + fRowCl->SetClass("AliTPCclusterMI"); + fRowCl->SetArray(1); + fRowCl->SetID(fParam->GetIndex(fSector, fRow)); + 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 + } // End of loop over sectors + + 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); - delete[] splitRows; - delete[] splitRowsRes; - Info("Digits2Clusters", "Number of found clusters : %d\n", nclusters); } -void AliTPCclustererMI::FindClusters() +void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC) { - //add virtual charge at the edge - for (Int_t i=0; i0){ @@ -720,45 +916,570 @@ void AliTPCclustererMI::FindClusters() } fBins[(fMaxPad+3)*fMaxTime+i] = amp0; } - -// memcpy(fResBins,fBins, fMaxBin*2); - memcpy(fResBins,fBins, fMaxBin); + memcpy(fResBins,fBins, fMaxBin*sizeof(Float_t)); + // + // // fNcluster=0; - //first loop - for "gold cluster" fLoop=1; Float_t *b=&fBins[-1]+2*fMaxTime; - Int_t crtime = Int_t((fParam->GetZLength()-AliTPCReconstructor::GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5); - + Int_t crtime = Int_t((fParam->GetZLength()-fRecoParam->GetCtgRange()*fRx)/fZWidth-fParam->GetNTBinsL1()-5); + Float_t minMaxCutAbs = fRecoParam->GetMinMaxCutAbs(); + Float_t minLeftRightCutAbs = fRecoParam->GetMinLeftRightCutAbs(); + Float_t minUpDownCutAbs = fRecoParam->GetMinUpDownCutAbs(); + Float_t minMaxCutSigma = fRecoParam->GetMinMaxCutSigma(); + Float_t minLeftRightCutSigma = fRecoParam->GetMinLeftRightCutSigma(); + Float_t minUpDownCutSigma = fRecoParam->GetMinUpDownCutSigma(); for (Int_t i=2*fMaxTime; iGetValue(fRow, i/fMaxTime); + // sigma cuts + if (b[0]GetFirstBin(); + // + UShort_t histo[kPedMax]; + memset(histo,0,kPedMax*sizeof(UShort_t)); + for (Int_t i=0; imax && i>firstBin) { + max = signal[i]; + maxPos = i; + } + if (signal[i]>kPedMax-1) continue; + histo[int(signal[i]+0.5)]++; + count0++; + } + // + for (Int_t i=1; ikPedMax) continue; + if (count06<0.6*count1){ + count06+=histo[median-idelta]; + mean06 +=histo[median-idelta]*(median-idelta); + rms06 +=histo[median-idelta]*(median-idelta)*(median-idelta); + count06+=histo[median+idelta]; + mean06 +=histo[median+idelta]*(median+idelta); + rms06 +=histo[median+idelta]*(median+idelta)*(median+idelta); + } + if (count09<0.9*count1){ + count09+=histo[median-idelta]; + mean09 +=histo[median-idelta]*(median-idelta); + rms09 +=histo[median-idelta]*(median-idelta)*(median-idelta); + count09+=histo[median+idelta]; + mean09 +=histo[median+idelta]*(median+idelta); + rms09 +=histo[median+idelta]*(median+idelta)*(median+idelta); + } + if (count10<0.95*count1){ + count10+=histo[median-idelta]; + mean +=histo[median-idelta]*(median-idelta); + rms +=histo[median-idelta]*(median-idelta)*(median-idelta); + count10+=histo[median+idelta]; + mean +=histo[median+idelta]*(median+idelta); + rms +=histo[median+idelta]*(median+idelta)*(median+idelta); + } + } + 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)); + rmsEvent = rms09; + // + pedestalEvent = median; + if (AliLog::GetDebugLevel("","AliTPCclustererMI")==0) return median; + // + UInt_t uid[3] = {UInt_t(id[0]),UInt_t(id[1]),UInt_t(id[2])}; + // + // Dump mean signal info + // + (*fDebugStreamer)<<"Signal"<< + "TimeStamp="<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]); + } + } + // + // + // + 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; i30.*TMath::Max(1.,rms06)){ + // + // + TGraph * graph =new TGraph(nchannels, dtime, dsignal); + // + // + // jumps left - right + Int_t njumps0=0; + Double_t deltaT0[2000]; + Double_t deltaA0[2000]; + Int_t lastJump0 = fRecoParam->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.,rms06) && + TMath::Abs(dsignal[itime]-dsignal[itime+1])>30.*TMath::Max(1.,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.,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.,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 (max-median>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="<fRecoParam->GetMaxNoise()) { + pedestalEvent+=1024.; + return 1024+median; // sign noisy channel in debug mode + } + return median; +} + + + +void AliTPCclustererMI::DumpHistos(){ + // + // Dump histogram information + // + if (!fAmplitudeHisto) return; + AliTPCROC * roc = AliTPCROC::Instance(); + 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 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)); + } } - */ } + + + +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; + } + + 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; +} +