X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCclustererMI.cxx;h=769730b547da36f6877c8625d4765d6c61d5f1d3;hb=1cb74071dd192ca26e97d2613cc05c9c8a91a398;hp=b0923dd14c8a6990194d8c7bb1abfb3bfc809b4d;hpb=679026e6d8485686f75155ed9a6e39ae3286522e;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCclustererMI.cxx b/TPC/AliTPCclustererMI.cxx index b0923dd14c8..769730b547d 100644 --- a/TPC/AliTPCclustererMI.cxx +++ b/TPC/AliTPCclustererMI.cxx @@ -25,6 +25,7 @@ // 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 @@ -43,6 +44,7 @@ #include #include #include +#include #include #include #include @@ -61,6 +63,7 @@ #include "AliTPCClustersRow.h" #include "AliTPCParam.h" #include "AliTPCRawStream.h" +#include "AliTPCRawStreamV3.h" #include "AliTPCRecoParam.h" #include "AliTPCReconstructor.h" #include "AliTPCcalibDB.h" @@ -79,7 +82,7 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar 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), @@ -95,13 +98,19 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar fInput(0), fOutput(0), fOutputArray(0), + fOutputClonesArray(0), fRowCl(0), fRowDig(0), fParam(0), fNcluster(0), + fNclusters(0), fDebugStreamer(0), fRecoParam(0), - fBDumpSignal(kFALSE) + fBDumpSignal(kFALSE), + fBClonesArray(kFALSE), + fAllBins(NULL), + fAllSigBins(NULL), + fAllNSigBins(NULL) { // // COSNTRUCTOR @@ -117,12 +126,33 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCParam* par, const AliTPCRecoPar fRecoParam = AliTPCReconstructor::GetRecoParam(); if (!fRecoParam) fRecoParam = AliTPCRecoParam::GetLowFluxParam(); } - fDebugStreamer = new TTreeSRedirector("TPCsignal.root"); + + if(AliTPCReconstructor::StreamLevel()>0) { + fDebugStreamer = new TTreeSRedirector("TPCsignal.root"); + } + // Int_t nPoints = fRecoParam->GetLastBin()-fRecoParam->GetFirstBin(); - fRowCl= new AliTPCClustersRow(); - fRowCl->SetClass("AliTPCclusterMI"); - fRowCl->SetArray(1); + 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) @@ -148,13 +178,19 @@ AliTPCclustererMI::AliTPCclustererMI(const AliTPCclustererMI ¶m) fInput(0), fOutput(0), fOutputArray(0), + fOutputClonesArray(0), fRowCl(0), fRowDig(0), fParam(0), fNcluster(0), + fNclusters(0), fDebugStreamer(0), fRecoParam(0), - fBDumpSignal(kFALSE) + fBDumpSignal(kFALSE), + fBClonesArray(kFALSE), + fAllBins(NULL), + fAllSigBins(NULL), + fAllNSigBins(NULL) { // // dummy @@ -180,6 +216,25 @@ AliTPCclustererMI::~AliTPCclustererMI(){ //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) @@ -203,10 +258,8 @@ void AliTPCclustererMI::SetOutput(TTree * tree) // if (!tree) return; fOutput= tree; - AliTPCClustersRow clrow; - AliTPCClustersRow *pclrow=&clrow; - clrow.SetClass("AliTPCclusterMI"); - clrow.SetArray(1); // to make Clones array + AliTPCClustersRow clrow, *pclrow=&clrow; + pclrow = new AliTPCClustersRow("AliTPCclusterMI"); fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200); } @@ -219,10 +272,10 @@ void AliTPCclustererMI::FillRow(){ // TObjArray // if (fOutput) fOutput->Fill(); - if (!fOutput){ + if (!fOutput && !fBClonesArray){ // if (!fOutputArray) fOutputArray = new TObjArray(fParam->GetNRowsTotal()); - if (fRowCl) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID()); + if (fRowCl && fRowCl->GetArray()->GetEntriesFast()>0) fOutputArray->AddAt(fRowCl->Clone(), fRowCl->GetID()); } } @@ -343,7 +396,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 @@ -393,8 +446,6 @@ AliTPCclusterMI &c) //unfolding 2 meani-=i0; meanj-=j0; - if (gDebug>4) - printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]); } @@ -431,7 +482,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]; @@ -522,8 +573,6 @@ void AliTPCclustererMI::UnfoldCluster(Float_t * matrix2[7], Float_t recmatrix[5] } } } - if (gDebug>4) - printf("%f\n", recmatrix[2][2]); } @@ -546,7 +595,7 @@ 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 rotated global coordinata @@ -579,10 +628,12 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t p // // // + AliTPCTransform *transform = AliTPCcalibDB::Instance()->GetTransform() ; if (!transform) { - AliFatal("Tranformations not in calibDB"); + AliFatal("Tranformations not in calibDB"); } + transform->SetCurrentRecoParam((AliTPCRecoParam*)fRecoParam); Double_t x[3]={c.GetRow(),c.GetPad(),c.GetTimeBin()}; Int_t i[1]={fSector}; transform->Transform(x,i,0,1); @@ -591,20 +642,23 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t p c.SetZ(x[2]); // // - if (!fRecoParam->GetBYMirror()){ - if (fSector%36>17){ - c.SetY(-c.GetY()); - } - } - if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) { c.SetType(-(c.GetType()+3)); //edge clusters } if (fLoop==2) c.SetType(100); if (!AcceptCluster(&c)) return; - TClonesArray * arr = fRowCl->GetArray(); - AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c); + // select output + TClonesArray * arr = 0; + AliTPCclusterMI * cl = 0; + + if(fBClonesArray==kFALSE) { + arr = fRowCl->GetArray(); + cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c); + } else { + cl = new ((*fOutputClonesArray)[fNclusters+fNcluster]) AliTPCclusterMI(c); + } + // if (fRecoParam->DumpSignal() &&matrix ) { // Int_t nbins=0; // Float_t *graph =0; @@ -620,8 +674,13 @@ void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c, Float_t * matrix, Int_t p } if (AliTPCReconstructor::StreamLevel()>1) { + Float_t xyz[3]; + cl->GetGlobalXYZ(xyz); (*fDebugStreamer)<<"Clusters"<< "Cl.="<5) { + AliInfo("Parameter Dumps"); + fParam->Dump(); + fRecoParam->Dump(); + } + AliTPCCalPad * gainTPC = AliTPCcalibDB::Instance()->GetPadGainFactor(); AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); AliSimDigits digarr, *dummy=&digarr; @@ -696,22 +764,161 @@ void AliTPCclustererMI::Digits2Clusters() Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3; Float_t gain = gainROC->GetValue(row,digarr.CurrentColumn()); Int_t bin = i*fMaxTime+j; - fBins[bin]=dig/gain; + if (gain>0){ + fBins[bin]=dig/gain; + }else{ + fBins[bin]=0; + } fSigBins[fNSigBins++]=bin; } while (digarr.Next()); digarr.ExpandTrackBuffer(); FindClusters(noiseROC); FillRow(); - fRowCl->GetArray()->Clear(); + fRowCl->GetArray()->Clear("C"); nclusters+=fNcluster; + delete[] fBins; delete[] fSigBins; } - + Info("Digits2Clusters", "Number of found clusters : %d", nclusters); } +void AliTPCclustererMI::ProcessSectorData(){ + // + // Process the data for the current sector + // + + AliTPCCalPad * pedestalTPC = AliTPCcalibDB::Instance()->GetPedestals(); + AliTPCCalPad * noiseTPC = AliTPCcalibDB::Instance()->GetPadNoise(); + AliTPCCalROC * pedestalROC = pedestalTPC->GetCalROC(fSector); // pedestal per given sector + AliTPCCalROC * noiseROC = noiseTPC->GetCalROC(fSector); // noise per given sector + //check the presence of the calibration + if (!noiseROC ||!pedestalROC ) { + AliError(Form("Missing calibration per sector\t%d\n",fSector)); + return; + } + Int_t nRows=fParam->GetNRow(fSector); + Bool_t calcPedestal = fRecoParam->GetCalcPedestal(); + Int_t zeroSup = fParam->GetZeroSup(); + // if (calcPedestal) { + if (kFALSE ) { + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad = fParam->GetNPads(fSector, iRow); + + for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) { + // + // Temporary fix for data production - !!!! MARIAN + // The noise calibration should take mean and RMS - currently the Gaussian fit used + // In case of double peak - the pad should be rejected + // + // Line mean - if more than given digits over threshold - make a noise calculation + // and pedestal substration + if (!calcPedestal && fAllBins[iRow][iPad*fMaxTime+0]<50) continue; + // + if (fAllBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data + Float_t *p = &fAllBins[iRow][iPad*fMaxTime+3]; + //Float_t pedestal = TMath::Median(fMaxTime, p); + Int_t id[3] = {fSector, iRow, iPad-3}; + // calib values + Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3); + Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3); + Double_t rmsEvent = rmsCalib; + Double_t pedestalEvent = pedestalCalib; + ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent); + if (rmsEventGetFirstBin()) + fAllBins[iRow][bin] = 0; + if (iTimeBin > fRecoParam->GetLastBin()) + fAllBins[iRow][bin] = 0; + if (fAllBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup) + fAllBins[iRow][bin] = 0; + if (fAllBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS + fAllBins[iRow][bin] = 0; + if (fAllBins[iRow][bin]) fAllSigBins[iRow][fAllNSigBins[iRow]++] = bin; + } + } + } + } + + if (AliTPCReconstructor::StreamLevel()>5) { + for (Int_t iRow = 0; iRow < nRows; iRow++) { + Int_t maxPad = fParam->GetNPads(fSector,iRow); + + for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) { + for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) { + Int_t bin = iPad*fMaxTime+iTimeBin; + Float_t signal = fAllBins[iRow][bin]; + if (AliTPCReconstructor::StreamLevel()>3 && signal>3) { + Double_t x[]={iRow,iPad-3,iTimeBin-3}; + Int_t i[]={fSector}; + AliTPCTransform trafo; + trafo.Transform(x,i,0,1); + Double_t gx[3]={x[0],x[1],x[2]}; + trafo.RotatedGlobal2Global(fSector,gx); + // fAllSigBins[iRow][fAllNSigBins[iRow]++] + Int_t rowsigBins = fAllNSigBins[iRow]; + Int_t first=fAllSigBins[iRow][0]; + Int_t last= 0; + // if (rowsigBins>0) fAllSigBins[iRow][fAllNSigBins[iRow]-1]; + + if (AliTPCReconstructor::StreamLevel()>5) { + (*fDebugStreamer)<<"Digits"<< + "sec="<SetID(fParam->GetIndex(fSector, fRow)); + if (fOutput) fOutput->GetBranch("Segment")->SetAddress(&fRowCl); + + fRx = fParam->GetPadRowRadii(fSector, fRow); + fPadLength = fParam->GetPadPitchLength(fSector, fRow); + fPadWidth = fParam->GetPadPitchWidth(); + fMaxPad = fParam->GetNPads(fSector,fRow); + fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after + + fBins = fAllBins[fRow]; + fSigBins = fAllSigBins[fRow]; + fNSigBins = fAllNSigBins[fRow]; + + FindClusters(noiseROC); + + FillRow(); + if(fBClonesArray == kFALSE) fRowCl->GetArray()->Clear("C"); + fNclusters += fNcluster; + + } // End of loop to find clusters +} + + void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) { //----------------------------------------------------------------- @@ -720,13 +927,211 @@ 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; + 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()); + } +} + + + + + +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(); AliTPCAltroMapping** mapping =AliTPCcalibDB::Instance()->GetMapping(); // AliTPCRawStream input(rawReader,(AliAltroMapping**)mapping); @@ -736,8 +1141,10 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) 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 = fRecoParam->GetLastBin() + 6; // add 3 virtual time bins before and 3 after const Int_t kNIS = fParam->GetNInnerSector(); @@ -746,37 +1153,23 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) fZWidth = fParam->GetZWidth(); Int_t zeroSup = fParam->GetZeroSup(); // - //alocate memory for sector - maximal case + // Clean-up // - Float_t** allBins = NULL; - Int_t** allSigBins = NULL; - Int_t* allNSigBins = 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]; - allSigBins = new Int_t*[nRowsMax]; - allNSigBins = new Int_t[nRowsMax]; for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { // Int_t maxBin = fMaxTime*(nPadsMax+6); // add 3 virtual pads before and 3 after - allBins[iRow] = new Float_t[maxBin]; - memset(allBins[iRow],0,sizeof(Float_t)*maxBin); - allSigBins[iRow] = new Int_t[maxBin]; - allNSigBins[iRow]=0; + 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 - //check the presence of the calibration - if (!noiseROC ||!pedestalROC ) { - AliError(Form("Missing calibration per sector\t%d\n",fSector)); - continue; - } Int_t nRows = 0; Int_t nDDLs = 0, indexDDL = 0; if (fSector < kNIS) { @@ -792,74 +1185,87 @@ 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 + input.Next(); + 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); - allNSigBins[iRow] = 0; + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + fAllNSigBins[iRow] = 0; } - // Loas the raw data for corresponding DDLs - rawReader->Reset(); - 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()) { 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){ - continue; + continue; } if (iRow < 0 || iRow >= nRows){ - AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !", + AliError(Form("Pad-row index (%d) outside the range (%d -> %d) !", iRow, 0, nRows -1)); - continue; + continue; } //pad Int_t iPad = input.GetPad(); if (iPad < 0 || iPad >= nPadsMax) { - AliError(Form("Pad index (%d) outside the range (%d -> %d) !", + AliError(Form("Pad index (%d) outside the range (%d -> %d) !", iPad, 0, nPadsMax-1)); - continue; + 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 < fRecoParam->GetFirstBin() || iTimeBin >= fRecoParam->GetLastBin()){ - continue; - AliFatal(Form("Timebin index (%d) outside the range (%d -> %d) !", + 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) { - Int_t bin = iPad*fMaxTime+iTimeBin; - allBins[iRow][bin] = signal/gain; - allSigBins[iRow][allNSigBins[iRow]++] = bin; + 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++; @@ -870,135 +1276,20 @@ void AliTPCclustererMI::Digits2Clusters(AliRawReader* rawReader) // // Now loop over rows and perform pedestal subtraction if (digCounter==0) continue; - // if (calcPedestal) { - if (kTRUE) { - for (Int_t iRow = 0; iRow < nRows; iRow++) { - Int_t maxPad; - if (fSector < kNIS) - maxPad = fParam->GetNPadsLow(iRow); - else - maxPad = fParam->GetNPadsUp(iRow); - - for (Int_t iPad = 3; iPad < maxPad + 3; iPad++) { - // - // Temporary fix for data production - !!!! MARIAN - // The noise calibration should take mean and RMS - currently the Gaussian fit used - // In case of double peak - the pad should be rejected - // - // Line mean - if more than given digits over threshold - make a noise calculation - // and pedestal substration - if (!calcPedestal && allBins[iRow][iPad*fMaxTime+0]<50) continue; - // - if (allBins[iRow][iPad*fMaxTime+0] <1 ) continue; // no data - Float_t *p = &allBins[iRow][iPad*fMaxTime+3]; - //Float_t pedestal = TMath::Median(fMaxTime, p); - Int_t id[3] = {fSector, iRow, iPad-3}; - // calib values - Double_t rmsCalib= noiseROC->GetValue(iRow,iPad-3); - Double_t pedestalCalib = pedestalROC->GetValue(iRow,iPad-3); - Double_t rmsEvent = rmsCalib; - Double_t pedestalEvent = pedestalCalib; - ProcesSignal(p, fMaxTime, id, rmsEvent, pedestalEvent); - if (rmsEventGetFirstBin()) - allBins[iRow][bin] = 0; - if (iTimeBin > AliTPCReconstructor::GetRecoParam()->GetLastBin()) - allBins[iRow][bin] = 0; - if (allBins[iRow][iPad*fMaxTime+iTimeBin] < zeroSup) - allBins[iRow][bin] = 0; - if (allBins[iRow][bin] < 3.0*rmsEvent) // 3 sigma cut on RMS - allBins[iRow][bin] = 0; - if (allBins[iRow][bin]) allSigBins[iRow][allNSigBins[iRow]++] = bin; - } - } - } - } - - if (AliTPCReconstructor::StreamLevel()>3) { - 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++) { - for (Int_t iTimeBin = 0; iTimeBin < fMaxTime; iTimeBin++) { - Int_t bin = iPad*fMaxTime+iTimeBin; - Float_t signal = allBins[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); - - (*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(); - if (fSector < kNIS) - fMaxPad = fParam->GetNPadsLow(fRow); - else - fMaxPad = fParam->GetNPadsUp(fRow); - fMaxBin = fMaxTime*(fMaxPad+6); // add 3 virtual pads before and 3 after - - fBins = allBins[fRow]; - fSigBins = allSigBins[fRow]; - fNSigBins = allNSigBins[fRow]; - - FindClusters(noiseROC); - FillRow(); - fRowCl->GetArray()->Clear(); - nclusters += fNcluster; - } // End of loop to find clusters + ProcessSectorData(); } // End of loop over sectors - - for (Int_t iRow = 0; iRow < nRowsMax; iRow++) { - delete [] allBins[iRow]; - delete [] allSigBins[iRow]; - } - delete [] allBins; - delete [] allSigBins; - delete [] allNSigBins; - + if (rawReader->GetEventId() && fOutput ){ - Info("Digits2Clusters", "File %s Event\t%d\tNumber of found clusters : %d\n", fOutput->GetName(),*(rawReader->GetEventId()), nclusters); - }else{ - Info("Digits2Clusters", "Event\t%d\tNumber of found clusters : %d\n",*(rawReader->GetEventId()), nclusters); - + 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(AliTPCCalROC * noiseROC) @@ -1023,6 +1314,7 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC) 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; @@ -1030,9 +1322,11 @@ void AliTPCclustererMI::FindClusters(AliTPCCalROC * noiseROC) //absolute custs if (b[0]GetFirstBin(); + Int_t firstBin = fRecoParam->GetFirstBin(); // UShort_t histo[kPedMax]; //memset(histo,0,kPedMax*sizeof(UShort_t)); @@ -1164,7 +1458,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="<kMin &&maxPos>AliTPCReconstructor::GetRecoParam()->GetFirstBin()) + if (AliTPCReconstructor::StreamLevel()>0) { + if (max-median>kMin &&maxPos>fRecoParam->GetFirstBin()) (*fDebugStreamer)<<"SignalB"<< // pads with signal "TimeStamp="<