X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=TPC%2FAliTPCdataQA.cxx;h=2cefe6129896e7f4fcc1ae3f7d79d3771009d561;hb=04b2edad5f2807838a2d8f737640a16e749fc992;hp=6b655bfcaa8399384af9af21c423b039d56c2dc7;hpb=336156cc1d6128122dd3cc39e1dd62d6aea6bc03;p=u%2Fmrichter%2FAliRoot.git diff --git a/TPC/AliTPCdataQA.cxx b/TPC/AliTPCdataQA.cxx index 6b655bfcaa8..2cefe612989 100644 --- a/TPC/AliTPCdataQA.cxx +++ b/TPC/AliTPCdataQA.cxx @@ -16,26 +16,90 @@ /* $Id$ */ +/* + June 2010 + + This update should solve two problems mainly: + * The vs event histograms have been limited to a fixed size for the + DQM. The 500k seemed to be a big size but is no longer so, so we + need to dynamically expand the range. The non-trivial point is that + we also have to do it for the copy owned by AliTPCQADataMakerRec. + * The amoreGui now remembers the x-range of the first visualization so + the trick of setting the relevant event range as the histogram is + filled no longer works. + + The fix is a bit crude but avoids creating a new histogram. Instead + the range is expanded (max events and events per bin is doubled) but + the number of bins is kept constant! In this way we can change just + double the max of the X-axis of the hist and rebin the data. The + same can easily be done for the copy owned by AliTPCQADataMakerRec. + + CAUTION: + If we change the number of bins we could crash the whole system + because ROOT does not create space for extra bins! (but we do not do + this). In that way it is a crude solution. + The rebinning in the code only works for an even number of bins. + + In addition to the above a bug in the reading of the config file was + also found and corrected. fAdcMax was set instead of fEventsPerBin. + + Finally cout was changes to AliInfo. + + February 2008 + + The code has been heavily modified so that now the RAW data is + "expanded" for each sector and stored in a big signal array. Then a + simple version of the code in AliTPCclustererMI is used to identify + the local maxima and these are then used for QA. This gives a better + estimate of the charge (both max and total) and also limits the + effect of noise. + + Implementation: + + In Update the RAW signals >= 3 ADC channels are stored in the arrays. + + There are 3 arrays: + Float_t** fAllBins 2d array [row][bin(pad, time)] ADC signal + Int_t** fAllSigBins 2d array [row][signal#] bin(with signal) + Int_t* fAllNSigBins; 1d array [row] Nsignals + + This is done sector by sector. + + When data from a new sector is encountered, the method + FindLocalMaxima is called on the data from the previous sector, and + the calibration/data objects are updated with the "cluster" + info. Finally the arrays are cleared. + + The requirements for a local maxima is: + Charge in bin is >= 5 ADC channels. + Charge in bin is larger than all the 8 neighboring bins. + At least one of the two pad neighbors has a signal. + At least one of the two time neighbors has a signal. -// stl includes -#include + Before accessing the data it is expected that the Analyse method is + called. This normalizes some of the data objects to per event or per + cluster. + If more data is passed to the class after Analyse has been called + the normalization is reversed and Analyse has to be called again. +*/ -using namespace std; //Root includes #include #include #include #include -#include -#include #include #include +#include +#include +#include //AliRoot includes #include "AliRawReader.h" #include "AliRawReaderRoot.h" #include "AliRawReaderDate.h" #include "AliTPCRawStream.h" +#include "AliTPCRawStreamV3.h" #include "AliTPCCalROC.h" #include "AliTPCROC.h" #include "AliMathBase.h" @@ -45,45 +109,98 @@ using namespace std; //date #include "event.h" #include "AliTPCCalPad.h" +#include "AliTPCPreprocessorOnline.h" //header file #include "AliTPCdataQA.h" +#include "AliLog.h" ClassImp(AliTPCdataQA) AliTPCdataQA::AliTPCdataQA() : /*FOLD00*/ - TH1F("TPCRAW","TPCRAW",100,0,100), fFirstTimeBin(60), fLastTimeBin(1000), fAdcMin(1), fAdcMax(100), - fOldRCUformat(kTRUE), - fROC(AliTPCROC::Instance()), fMapping(NULL), fPedestal(0), fNoise(0), + fNLocalMaxima(0), fMaxCharge(0), fMeanCharge(0), fNoThreshold(0), - fOverThreshold0(0), - fOverThreshold5(0), + fNTimeBins(0), + fNPads(0), + fTimePosition(0), fOverThreshold10(0), fOverThreshold20(0), fOverThreshold30(0), - fEventCounter(0) + fHistQVsTimeSideA(0), + fHistQVsTimeSideC(0), + fHistQMaxVsTimeSideA(0), + fHistQMaxVsTimeSideC(0), + fHistOccupancyVsEvent(0), + fHistNclustersVsEvent(0), + fEventCounter(0), + fIsAnalysed(kFALSE), + fMaxEvents(500000), // Max events for event histograms + fEventsPerBin(1000), // Events per bin for event histograms + fSignalCounter(0), // Signal counter + fClusterCounter(0), // Cluster counter + fAllBins(0), + fAllSigBins(0), + fAllNSigBins(0), + fRowsMax(0), + fPadsMax(0), + fTimeBinsMax(0) { // // default constructor // - - fSectorLast = -1; - fRowLast = 0; - fPadLast = 0; - fTimeBinLast = 0; - fSignalLast = 0; - fNAboveThreshold = 0; } +//_____________________________________________________________________ +AliTPCdataQA::AliTPCdataQA(AliRecoParam::EventSpecie_t es) : +fFirstTimeBin(60), +fLastTimeBin(1000), +fAdcMin(1), +fAdcMax(100), +fMapping(NULL), +fPedestal(0), +fNoise(0), +fNLocalMaxima(0), +fMaxCharge(0), +fMeanCharge(0), +fNoThreshold(0), +fNTimeBins(0), +fNPads(0), +fTimePosition(0), +fOverThreshold10(0), +fOverThreshold20(0), +fOverThreshold30(0), +fHistQVsTimeSideA(0), +fHistQVsTimeSideC(0), +fHistQMaxVsTimeSideA(0), +fHistQMaxVsTimeSideC(0), +fHistOccupancyVsEvent(0), +fHistNclustersVsEvent(0), +fEventCounter(0), +fIsAnalysed(kFALSE), +fMaxEvents(500000), +fEventsPerBin(1000), +fSignalCounter(0), +fClusterCounter(0), +fAllBins(0), +fAllSigBins(0), +fAllNSigBins(0), +fRowsMax(0), +fPadsMax(0), +fTimeBinsMax(0) +{ +// ctor creating the histogram + char * name = Form("TPCRAW_%s", AliRecoParam::GetEventSpecieName(es)) ; + TH1F(name, name,100,0,100) ; +} //_____________________________________________________________________ AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/ @@ -92,16 +209,136 @@ AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/ fLastTimeBin(ped.GetLastTimeBin()), fAdcMin(ped.GetAdcMin()), fAdcMax(ped.GetAdcMax()), - fOldRCUformat(ped.fOldRCUformat), - fROC(AliTPCROC::Instance()), - fMapping(NULL) + fMapping(NULL), + fPedestal(0), + fNoise(0), + fNLocalMaxima(0), + fMaxCharge(0), + fMeanCharge(0), + fNoThreshold(0), + fNTimeBins(0), + fNPads(0), + fTimePosition(0), + fOverThreshold10(0), + fOverThreshold20(0), + fOverThreshold30(0), + fHistQVsTimeSideA(0), + fHistQVsTimeSideC(0), + fHistQMaxVsTimeSideA(0), + fHistQMaxVsTimeSideC(0), + fHistOccupancyVsEvent(0), + fHistNclustersVsEvent(0), + fEventCounter(ped.GetEventCounter()), + fIsAnalysed(ped.GetIsAnalysed()), + fMaxEvents(ped.GetMaxEvents()), + fEventsPerBin(ped.GetEventsPerBin()), + fSignalCounter(ped.GetSignalCounter()), + fClusterCounter(ped.GetClusterCounter()), + fAllBins(0), + fAllSigBins(0), + fAllNSigBins(0), + fRowsMax(0), + fPadsMax(0), + fTimeBinsMax(0) { // // copy constructor // - + if(ped.GetNLocalMaxima()) + fNLocalMaxima = new AliTPCCalPad(*ped.GetNLocalMaxima()); + if(ped.GetMaxCharge()) + fMaxCharge = new AliTPCCalPad(*ped.GetMaxCharge()); + if(ped.GetMeanCharge()) + fMeanCharge = new AliTPCCalPad(*ped.GetMeanCharge()); + if(ped.GetNoThreshold()) + fNoThreshold = new AliTPCCalPad(*ped.GetNoThreshold()); + if(ped.GetNTimeBins()) + fNTimeBins = new AliTPCCalPad(*ped.GetNTimeBins()); + if(ped.GetNPads()) + fNPads = new AliTPCCalPad(*ped.GetNPads()); + if(ped.GetTimePosition()) + fTimePosition = new AliTPCCalPad(*ped.GetTimePosition()); + if(ped.GetOverThreshold10()) + fOverThreshold10 = new AliTPCCalPad(*ped.GetOverThreshold10()); + if(ped.GetOverThreshold20()) + fOverThreshold20 = new AliTPCCalPad(*ped.GetOverThreshold20()); + if(ped.GetOverThreshold30()) + fOverThreshold30 = new AliTPCCalPad(*ped.GetOverThreshold30()); + if(ped.GetHistQVsTimeSideA()) { + fHistQVsTimeSideA = new TProfile(*ped.GetHistQVsTimeSideA()); + fHistQVsTimeSideA->SetDirectory(0); + } + if(ped.GetHistQVsTimeSideC()) { + fHistQVsTimeSideC = new TProfile(*ped.GetHistQVsTimeSideC()); + fHistQVsTimeSideC->SetDirectory(0); + } + if(ped.GetHistQMaxVsTimeSideA()) { + fHistQMaxVsTimeSideA = new TProfile(*ped.GetHistQMaxVsTimeSideA()); + fHistQMaxVsTimeSideA->SetDirectory(0); + } + if(ped.GetHistQMaxVsTimeSideC()) { + fHistQMaxVsTimeSideC = new TProfile(*ped.GetHistQMaxVsTimeSideC()); + fHistQMaxVsTimeSideC->SetDirectory(0); + } + if(ped.GetHistOccupancyVsEventConst()) { + fHistOccupancyVsEvent = new TH1F(*ped.GetHistOccupancyVsEventConst()); + fHistOccupancyVsEvent->SetDirectory(0); + } + if(ped.GetHistNclustersVsEventConst()) { + fHistNclustersVsEvent = new TH1F(*ped.GetHistNclustersVsEventConst()); + fHistNclustersVsEvent->SetDirectory(0); + } } +//_____________________________________________________________________ +AliTPCdataQA::AliTPCdataQA(const TMap *config) : /*FOLD00*/ + TH1F("TPCRAW","TPCRAW",100,0,100), + fFirstTimeBin(60), + fLastTimeBin(1000), + fAdcMin(1), + fAdcMax(100), + fMapping(NULL), + fPedestal(0), + fNoise(0), + fNLocalMaxima(0), + fMaxCharge(0), + fMeanCharge(0), + fNoThreshold(0), + fNTimeBins(0), + fNPads(0), + fTimePosition(0), + fOverThreshold10(0), + fOverThreshold20(0), + fOverThreshold30(0), + fHistQVsTimeSideA(0), + fHistQVsTimeSideC(0), + fHistQMaxVsTimeSideA(0), + fHistQMaxVsTimeSideC(0), + fHistOccupancyVsEvent(0), + fHistNclustersVsEvent(0), + fEventCounter(0), + fIsAnalysed(kFALSE), + fMaxEvents(500000), + fEventsPerBin(1000), + fSignalCounter(0), + fClusterCounter(0), + fAllBins(0), + fAllSigBins(0), + fAllNSigBins(0), + fRowsMax(0), + fPadsMax(0), + fTimeBinsMax(0) +{ + // + // default constructor + // + if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi(); + if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi(); + if (config->GetValue("AdcMin")) fAdcMin = ((TObjString*)config->GetValue("AdcMin"))->GetString().Atoi(); + if (config->GetValue("AdcMax")) fAdcMax = ((TObjString*)config->GetValue("AdcMax"))->GetString().Atoi(); + if (config->GetValue("MaxEvents")) fMaxEvents = ((TObjString*)config->GetValue("MaxEvents"))->GetString().Atoi(); + if (config->GetValue("EventsPerBin")) fEventsPerBin = ((TObjString*)config->GetValue("EventsPerBin"))->GetString().Atoi(); +} //_____________________________________________________________________ AliTPCdataQA& AliTPCdataQA::operator = (const AliTPCdataQA &source) @@ -124,68 +361,272 @@ AliTPCdataQA::~AliTPCdataQA() /*FOLD00*/ // // do not delete fMapping, because we do not own it. + // do not delete fMapping, because we do not own it. + // do not delete fNoise and fPedestal, because we do not own them. + + delete fNLocalMaxima; + delete fMaxCharge; + delete fMeanCharge; + delete fNoThreshold; + delete fNTimeBins; + delete fNPads; + delete fTimePosition; + delete fOverThreshold10; + delete fOverThreshold20; + delete fOverThreshold30; + delete fHistQVsTimeSideA; + delete fHistQVsTimeSideC; + delete fHistQMaxVsTimeSideA; + delete fHistQMaxVsTimeSideC; + delete fHistOccupancyVsEvent; + delete fHistNclustersVsEvent; + + for (Int_t iRow = 0; iRow < fRowsMax; iRow++) { + delete [] fAllBins[iRow]; + delete [] fAllSigBins[iRow]; + } + delete [] fAllBins; + delete [] fAllSigBins; + delete [] fAllNSigBins; +} + +//_____________________________________________________________________ +TH1F* AliTPCdataQA::GetHistOccupancyVsEvent() +{ + // + // Create Occupancy vs event histogram + // (we create this histogram differently then the other histograms + // because this we want to be able to access and copy + // from AliTPCQAMakerRec before it normally would be created) + // + if(!fHistOccupancyVsEvent) { + Int_t nBins = fMaxEvents/fEventsPerBin; + fHistOccupancyVsEvent = new TH1F("hOccupancyVsEvent", "Occupancy vs event number (~time); Event number; Occupancy", nBins, 0, nBins*fEventsPerBin); + fHistOccupancyVsEvent->SetDirectory(0); + } + + return fHistOccupancyVsEvent; } +//_____________________________________________________________________ +TH1F* AliTPCdataQA::GetHistNclustersVsEvent() +{ + // + // Create Nclusters vs event histogram + // (we create this histogram differently then the other histograms + // because this we want to be able to access and copy + // from AliTPCQAMakerRec before it normally would be created) + // + if(!fHistNclustersVsEvent) { + Int_t nBins = fMaxEvents/fEventsPerBin; + fHistNclustersVsEvent = new TH1F("hNclustersVsEvent", "Nclusters vs event number (~time); Event number; Nclusters per event", nBins, 0, nBins*fEventsPerBin); + fHistNclustersVsEvent->SetDirectory(0); + } + + return fHistNclustersVsEvent; +} +//_____________________________________________________________________ +void AliTPCdataQA::UpdateEventHistograms() +{ + // Update histograms that display occupancy and + // number of clusters as a function of number of + // events + if (!fHistOccupancyVsEvent) + GetHistOccupancyVsEvent(); + if (!fHistNclustersVsEvent) + GetHistNclustersVsEvent(); + + if(fEventCounter > fMaxEvents) { + + // we have to expand the histogram to handle the larger number of + // events. The way it is done now is to double the range and the + // number of events per bin (so the number of histogram bins stays + // constant) + fEventsPerBin *= 2; + fMaxEvents *= 2; + + // Change histogram limits + const Int_t nBins = fHistOccupancyVsEvent->GetXaxis()->GetNbins(); + fHistOccupancyVsEvent->GetXaxis()->Set(nBins, fHistOccupancyVsEvent->GetXaxis()->GetNbins(), fMaxEvents); + fHistNclustersVsEvent->GetXaxis()->Set(nBins, fHistNclustersVsEvent->GetXaxis()->GetNbins(), fMaxEvents); + + // Rebin the histogram + for(Int_t bin = 1; bin <= nBins; bin+=2) { + + Int_t newBin = TMath::Nint(Float_t(bin+1)/2.0); + Float_t newContent = (fHistOccupancyVsEvent->GetBinContent(bin)+ + fHistOccupancyVsEvent->GetBinContent(bin+1))/2.0; + fHistOccupancyVsEvent->SetBinContent(newBin, newContent); + + newContent = (fHistNclustersVsEvent->GetBinContent(bin)+ + fHistNclustersVsEvent->GetBinContent(bin+1))/2.0; + fHistNclustersVsEvent->SetBinContent(newBin, newContent); + } + + // Set the remaining bins to 0 + Int_t lastHalf = nBins/2 +1; + for(Int_t bin = lastHalf; bin <= nBins; bin++) { + + fHistOccupancyVsEvent->SetBinContent(bin, 0); + fHistNclustersVsEvent->SetBinContent(bin, 0); + } + + // In this case we should nut update but wait untill the new + // number of events per bin is reached! + return; + } + + const Int_t bin = TMath::Nint(Float_t(fEventCounter)/fEventsPerBin); + + Float_t averageOccupancy = + Float_t(fSignalCounter)/fEventsPerBin/(fLastTimeBin - fFirstTimeBin +1.0) + / 570132.0; // 570,132 is number of pads + fHistOccupancyVsEvent->SetBinContent(bin, averageOccupancy); + fSignalCounter = 0; + + Float_t averageNclusters = + Float_t(fClusterCounter)/fEventsPerBin; + fHistNclustersVsEvent->SetBinContent(bin, averageNclusters); + fClusterCounter = 0; +} + +//_____________________________________________________________________ +Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStreamV3 *const rawStreamV3) +{ + // + // Event Processing loop - AliTPCRawStreamV3 + // + Bool_t withInput = kFALSE; + Int_t nSignals = 0; + Int_t lastSector = -1; + + while ( rawStreamV3->NextDDL() ){ + + while ( rawStreamV3->NextChannel() ){ + + Int_t iSector = rawStreamV3->GetSector(); // current sector + Int_t iRow = rawStreamV3->GetRow(); // current row + Int_t iPad = rawStreamV3->GetPad(); // current pad + if (iRow<0 || iPad<0) continue; + // Call local maxima finder if the data is in a new sector + if(iSector != lastSector) { + + if(nSignals>0) + FindLocalMaxima(lastSector); + + CleanArrays(); + lastSector = iSector; + nSignals = 0; + } + + while ( rawStreamV3->NextBunch() ){ + + Int_t startTbin = (Int_t)rawStreamV3->GetStartTimeBin(); + Int_t bunchlength = (Int_t)rawStreamV3->GetBunchLength(); + const UShort_t *sig = rawStreamV3->GetSignals(); + + for (Int_t iTimeBin = 0; iTimeBin=0&&nSignals>0) + FindLocalMaxima(lastSector); + + return withInput; +} //_____________________________________________________________________ -Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast) +Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *const rawReader) +{ + // + // Event processing loop - AliRawReader + // + AliTPCRawStreamV3 rawStreamV3(rawReader, (AliAltroMapping**)fMapping); + Bool_t res=ProcessEvent(&rawStreamV3); + if(res) { + fEventCounter++; // only increment event counter if there is TPC data + + if(fEventCounter%fEventsPerBin==0) + UpdateEventHistograms(); + } + return res; +} + +//_____________________________________________________________________ +Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *const rawStreamFast) { // // Event Processing loop - AliTPCRawStream // Bool_t withInput = kFALSE; + Int_t nSignals = 0; + Int_t lastSector = -1; while ( rawStreamFast->NextDDL() ){ - while ( rawStreamFast->NextChannel() ){ - Int_t isector = rawStreamFast->GetSector(); // current sector - Int_t iRow = rawStreamFast->GetRow(); // current row - Int_t iPad = rawStreamFast->GetPad(); // current pad - - while ( rawStreamFast->NextBunch() ){ - Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin(); - Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin(); - - for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){ - Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin]; - Update(isector,iRow,iPad,iTimeBin+1,signal); - withInput = kTRUE; - } - } + while ( rawStreamFast->NextChannel() ){ + + Int_t iSector = rawStreamFast->GetSector(); // current sector + Int_t iRow = rawStreamFast->GetRow(); // current row + Int_t iPad = rawStreamFast->GetPad(); // current pad + // Call local maxima finder if the data is in a new sector + if(iSector != lastSector) { + + if(nSignals>0) + FindLocalMaxima(lastSector); + + CleanArrays(); + lastSector = iSector; + nSignals = 0; + } + + while ( rawStreamFast->NextBunch() ){ + Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin(); + Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin(); + + for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){ + Float_t signal = rawStreamFast->GetSignals()[iTimeBin-startTbin]; + nSignals += Update(iSector,iRow,iPad,iTimeBin+1,signal); + withInput = kTRUE; + } } + } } - + return withInput; } //_____________________________________________________________________ -Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *rawReader) +Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *const rawReader) { // // Event processing loop - AliRawReader // - fSectorLast = -1; - AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping); - Bool_t res=ProcessEventFast(rawStreamFast); + AliTPCRawStreamFast rawStreamFast(rawReader, (AliAltroMapping**)fMapping); + Bool_t res=ProcessEventFast(&rawStreamFast); if(res) fEventCounter++; // only increment event counter if there is TPC data // otherwise Analyse (called in QA) fails - delete rawStreamFast; return res; } //_____________________________________________________________________ -Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream) +Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *const rawStream) { // // Event Processing loop - AliTPCRawStream // - rawStream->SetOldRCUFormat(fOldRCUformat); - Bool_t withInput = kFALSE; + Int_t nSignals = 0; + Int_t lastSector = -1; while (rawStream->Next()) { @@ -195,23 +636,39 @@ Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream) Int_t iTimeBin = rawStream->GetTime(); // current time bin Float_t signal = rawStream->GetSignal(); // current ADC signal - Update(iSector,iRow,iPad,iTimeBin,signal); - withInput = kTRUE; + // Call local maxima finder if the data is in a new sector + if(iSector != lastSector) { + + if(nSignals>0) + FindLocalMaxima(lastSector); + + CleanArrays(); + lastSector = iSector; + nSignals = 0; + } + + // Sometimes iRow==-1 if there is problems to read the data + if(iRow>=0) { + nSignals += Update(iSector,iRow,iPad,iTimeBin,signal); + withInput = kTRUE; + } } + if (lastSector>=0&&nSignals>0) + FindLocalMaxima(lastSector); + return withInput; } //_____________________________________________________________________ -Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader) +Bool_t AliTPCdataQA::ProcessEventOld(AliRawReader *const rawReader) { // // Event processing loop - AliRawReader // // if fMapping is NULL the rawstream will crate its own mapping - fSectorLast = -1; AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping); rawReader->Select("TPC"); Bool_t res = ProcessEvent(&rawStream); @@ -224,15 +681,14 @@ Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader) //_____________________________________________________________________ -Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *event) +Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *const event) { // // process date event // - AliRawReader *rawReader = new AliRawReaderDate((void*)event); - Bool_t result=ProcessEvent(rawReader); - delete rawReader; + AliRawReaderDate rawReader((void*)event); + Bool_t result=ProcessEvent(&rawReader); return result; } @@ -268,135 +724,260 @@ void AliTPCdataQA::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t //_____________________________________________________________________ -Int_t AliTPCdataQA::Update(const Int_t icsector, /*FOLD00*/ - const Int_t icRow, - const Int_t icPad, - const Int_t icTimeBin, - const Float_t csignal) +Int_t AliTPCdataQA::Update(const Int_t iSector, /*FOLD00*/ + const Int_t iRow, + const Int_t iPad, + const Int_t iTimeBin, + Float_t signal) { // // Signal filling method // - if (icTimeBinfLastTimeBin) return 0; - + + // + // Define the calibration objects the first time Update is called + // NB! This has to be done first even if the data is rejected by the time + // cut to make sure that the objects are available in Analyse + // + if (!fNLocalMaxima) fNLocalMaxima = new AliTPCCalPad("NLocalMaxima","NLocalMaxima"); if (!fMaxCharge) fMaxCharge = new AliTPCCalPad("MaxCharge","MaxCharge"); if (!fMeanCharge) fMeanCharge = new AliTPCCalPad("MeanCharge","MeanCharge"); if (!fNoThreshold) fNoThreshold = new AliTPCCalPad("NoThreshold","NoThreshold"); - if (!fOverThreshold0) fOverThreshold0 = new AliTPCCalPad("OverThreshold0","OverThreshold0"); - if (!fOverThreshold5) fOverThreshold5 = new AliTPCCalPad("OverThreshold5","OverThreshold5"); + if (!fNTimeBins) fNTimeBins = new AliTPCCalPad("NTimeBins","NTimeBins"); + if (!fNPads) fNPads = new AliTPCCalPad("NPads","NPads"); + if (!fTimePosition) fTimePosition = new AliTPCCalPad("TimePosition","TimePosition"); if (!fOverThreshold10) fOverThreshold10 = new AliTPCCalPad("OverThreshold10","OverThreshold10"); if (!fOverThreshold20) fOverThreshold20 = new AliTPCCalPad("OverThreshold20","OverThreshold20"); if (!fOverThreshold30) fOverThreshold30 = new AliTPCCalPad("OverThreshold30","OverThreshold30"); - // + if (!fHistQVsTimeSideA) { + fHistQVsTimeSideA = new TProfile("hQVsTimeSideA", "Q vs time (side A); Time [Timebin]; Q [ADC ch]", 100, 0, 1000); + fHistQVsTimeSideA->SetDirectory(0); + } + if (!fHistQVsTimeSideC) { + fHistQVsTimeSideC = new TProfile("hQVsTimeSideC", "Q vs time (side C); Time [Timebin]; Q [ADC ch]", 100, 0, 1000); + fHistQVsTimeSideC->SetDirectory(0); + } + if (!fHistQMaxVsTimeSideA) { + fHistQMaxVsTimeSideA = new TProfile("hQMaxVsTimeSideA", "Q_{MAX} vs time (side A); Time [Timebin]; Q_{MAX} [ADC ch]", 100, 0, 1000); + fHistQMaxVsTimeSideA->SetDirectory(0); + } + if (!fHistQMaxVsTimeSideC) { + fHistQMaxVsTimeSideC = new TProfile("hQMaxVsTimeSideC", "Q_{MAX} vs time (side C); Time [Timebin]; Q_{MAX} [ADC ch]", 100, 0, 1000); + fHistQMaxVsTimeSideC->SetDirectory(0); + } - Int_t signal = Int_t(csignal); + // Make the arrays for expanding the data - // if pedestal calibrations are loaded subtract pedestals - if(fPedestal) { + if (!fAllBins) + MakeArrays(); - Int_t pedestal = Int_t(fPedestal->GetCalROC(icsector)->GetValue(icRow, icPad)); - if(pedestal<10 || pedestal>90) - return 0; - signal -= pedestal; + // + // If Analyse has been previously called we need now to denormalize the data + // as more data is coming + // + if(fIsAnalysed == kTRUE) { + + const Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1; + const Float_t denormalization = Float_t(fEventCounter * nTimeBins); + fNoThreshold->Multiply(denormalization); + + fMeanCharge->Multiply(fNLocalMaxima); + fMaxCharge->Multiply(fNLocalMaxima); + fNTimeBins->Multiply(fNLocalMaxima); + fNPads->Multiply(fNLocalMaxima); + fTimePosition->Multiply(fNLocalMaxima); + fIsAnalysed = kFALSE; } + // + // TimeBin cut + // + if (iTimeBinfLastTimeBin) return 0; + + // if pedestal calibrations are loaded subtract pedestals + if(fPedestal) { - if (signal >= 0) { - - Float_t count = fNoThreshold->GetCalROC(icsector)->GetValue(icRow, icPad); - fNoThreshold->GetCalROC(icsector)->SetValue(icRow, icPad,count+1); + Float_t ped = fPedestal->GetCalROC(iSector)->GetValue(iRow, iPad); + // Don't use data from pads where pedestals are abnormally small or large + if(ped<10 || ped>90) + return 0; + signal -= ped; } - + + // In fNoThreshold we fill all data to estimate the ZS volume + Float_t count = fNoThreshold->GetCalROC(iSector)->GetValue(iRow, iPad); + fNoThreshold->GetCalROC(iSector)->SetValue(iRow, iPad,count+1); + // Require at least 3 ADC channels - if (signal < 3) + if (signal < 3.0) return 0; // if noise calibrations are loaded require at least 3*sigmaNoise if(fNoise) { - - Float_t noise = fNoise->GetCalROC(icsector)->GetValue(icRow, icPad); - - if(signalGetCalROC(iSector)->GetValue(iRow, iPad); + + if(signal < noise*3.0) return 0; } + // - // this signal is ok - now see if the previous signal was connected - // this is a bit ugly since the standard decoder goes down in time bins - // (10, 9, 8..) while the fast HLT decoder goes up in time bins (1, 2, 3..) + // This signal is ok and we store it in the cluster map // - if(fSectorLast==icsector && fRowLast==icRow && fPadLast==icPad && - fTimeBinLast==icTimeBin+1 || fTimeBinLast==icTimeBin-1) - fNAboveThreshold++; - else - fNAboveThreshold = 1; - - if(fNAboveThreshold==2) { - - // - // This is the only special case, because before we did not know if we - // should store the information - // - UpdateSignalHistograms(fSectorLast, fRowLast, fPadLast, fTimeBinLast, - fSignalLast); - } - - // keep the information for the next signal - fSectorLast = icsector; - fRowLast = icRow; - fPadLast = icPad; - fTimeBinLast = icTimeBin; - fSignalLast = signal; - - if(fNAboveThreshold==1) // we don't know if this information should be stored - return 1; - - UpdateSignalHistograms(fSectorLast, fRowLast, fPadLast, fTimeBinLast, - fSignalLast); - return 1; + SetExpandDigit(iRow, iPad, iTimeBin, signal); + + fSignalCounter++; + + return 1; // signal was accepted } + //_____________________________________________________________________ -void AliTPCdataQA::UpdateSignalHistograms(const Int_t icsector, /*FOLD00*/ - const Int_t icRow, - const Int_t icPad, - const Int_t icTimeBin, - const Float_t signal) +void AliTPCdataQA::FindLocalMaxima(const Int_t iSector) { // - // Signal filling method + // This method is called after the data from each sector has been + // exapanded into an array + // Loop over the signals and identify local maxima and fill the + // calibration objects with the information // + + Int_t nLocalMaxima = 0; + const Int_t maxTimeBin = fTimeBinsMax+4; // Used to step between neighboring pads + // Because we have tha pad-time data in a + // 1d array + + for (Int_t iRow = 0; iRow < fRowsMax; iRow++) { + + Float_t* allBins = fAllBins[iRow]; + Int_t* sigBins = fAllSigBins[iRow]; + const Int_t nSigBins = fAllNSigBins[iRow]; + + for (Int_t iSig = 0; iSig < nSigBins; iSig++) { + + Int_t bin = sigBins[iSig]; + Float_t *b = &allBins[bin]; + + // + // Now we check if this is a local maximum + // + + Float_t qMax = b[0]; + + // First check that the charge is bigger than the threshold + if (qMax<5) + continue; + + // Require at least one neighboring pad with signal + if (b[-maxTimeBin]+b[maxTimeBin]<=0) continue; + + // Require at least one neighboring time bin with signal + if (b[-1]+b[1]<=0) continue; + + // + // Check that this is a local maximum + // Note that the checking is done so that if 2 charges has the same + // qMax then only 1 cluster is generated + // (that is why there is BOTH > and >=) + // + if (b[-maxTimeBin] >= qMax) continue; + if (b[-1 ] >= qMax) continue; + if (b[+maxTimeBin] > qMax) continue; + if (b[+1 ] > qMax) continue; + if (b[-maxTimeBin-1] >= qMax) continue; + if (b[+maxTimeBin-1] >= qMax) continue; + if (b[+maxTimeBin+1] > qMax) continue; + if (b[-maxTimeBin+1] >= qMax) continue; + + // + // Now we accept the local maximum and fill the calibration/data objects + // + nLocalMaxima++; + + Int_t iPad, iTimeBin; + GetPadAndTimeBin(bin, iPad, iTimeBin); + + Float_t count = fNLocalMaxima->GetCalROC(iSector)->GetValue(iRow, iPad); + fNLocalMaxima->GetCalROC(iSector)->SetValue(iRow, iPad, count+1); + + count = fTimePosition->GetCalROC(iSector)->GetValue(iRow, iPad); + fTimePosition->GetCalROC(iSector)->SetValue(iRow, iPad, count+iTimeBin); + + Float_t charge = fMaxCharge->GetCalROC(iSector)->GetValue(iRow, iPad); + fMaxCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qMax); + + if(qMax>=10) { + count = fOverThreshold10->GetCalROC(iSector)->GetValue(iRow, iPad); + fOverThreshold10->GetCalROC(iSector)->SetValue(iRow, iPad, count+1); + } + if(qMax>=20) { + count = fOverThreshold20->GetCalROC(iSector)->GetValue(iRow, iPad); + fOverThreshold20->GetCalROC(iSector)->SetValue(iRow, iPad, count+1); + } + if(qMax>=30) { + count = fOverThreshold30->GetCalROC(iSector)->GetValue(iRow, iPad); + fOverThreshold30->GetCalROC(iSector)->SetValue(iRow, iPad, count+1); + } + + // + // Calculate the total charge as the sum over the region: + // + // o o o o o + // o i i i o + // o i C i o + // o i i i o + // o o o o o + // + // with qmax at the center C. + // + // The inner charge (i) we always add, but we only add the outer + // charge (o) if the neighboring inner bin (i) has a signal. + // + Int_t minP = 0, maxP = 0, minT = 0, maxT = 0; + Float_t qTot = qMax; + for(Int_t i = -1; i<=1; i++) { + for(Int_t j = -1; j<=1; j++) { + + if(i==0 && j==0) + continue; + + Float_t charge1 = GetQ(b, i, j, maxTimeBin, minT, maxT, minP, maxP); + qTot += charge1; + if(charge1>0) { + // see if the next neighbor is also above threshold + if(i*j==0) { + qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP); + } else { + // we are in a diagonal corner + qTot += GetQ(b, i, 2*j, maxTimeBin, minT, maxT, minP, maxP); + qTot += GetQ(b, 2*i, j, maxTimeBin, minT, maxT, minP, maxP); + qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP); + } + } + } + } + + charge = fMeanCharge->GetCalROC(iSector)->GetValue(iRow, iPad); + fMeanCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qTot); + + count = fNTimeBins->GetCalROC(iSector)->GetValue(iRow, iPad); + fNTimeBins->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxT-minT+1); + + count = fNPads->GetCalROC(iSector)->GetValue(iRow, iPad); + fNPads->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxP-minP+1); + + if((iSector%36)<18) { // A side + fHistQVsTimeSideA->Fill(iTimeBin, qTot); + fHistQMaxVsTimeSideA->Fill(iTimeBin, qMax); + } else { + fHistQVsTimeSideC->Fill(iTimeBin, qTot); + fHistQMaxVsTimeSideC->Fill(iTimeBin, qMax); + } + } // end loop over signals + } // end loop over rows - { - Float_t charge = fMeanCharge->GetCalROC(icsector)->GetValue(icRow, icPad); - fMeanCharge->GetCalROC(icsector)->SetValue(icRow, icPad, charge + signal); - } - - if (signal>fMaxCharge->GetCalROC(icsector)->GetValue(icRow, icPad)){ - fMaxCharge->GetCalROC(icsector)->SetValue(icRow, icPad,signal); - } - - if (signal>0){ - Float_t count = fOverThreshold0->GetCalROC(icsector)->GetValue(icRow, icPad); - fOverThreshold0->GetCalROC(icsector)->SetValue(icRow, icPad,count+1); - }; - // - if (signal>5){ - Float_t count = fOverThreshold5->GetCalROC(icsector)->GetValue(icRow, icPad); - fOverThreshold5->GetCalROC(icsector)->SetValue(icRow, icPad,count+1); - }; - if (signal>10){ - Float_t count = fOverThreshold10->GetCalROC(icsector)->GetValue(icRow, icPad); - fOverThreshold10->GetCalROC(icsector)->SetValue(icRow, icPad,count+1); - }; - if (signal>20){ - Float_t count = fOverThreshold20->GetCalROC(icsector)->GetValue(icRow, icPad); - fOverThreshold20->GetCalROC(icsector)->SetValue(icRow, icPad,count+1); - }; - if (signal>30){ - Float_t count = fOverThreshold30->GetCalROC(icsector)->GetValue(icRow, icPad); - fOverThreshold30->GetCalROC(icsector)->SetValue(icRow, icPad,count+1); - }; + fClusterCounter += nLocalMaxima; } //_____________________________________________________________________ @@ -406,26 +987,214 @@ void AliTPCdataQA::Analyse() // Calculate calibration constants // - cout << "Analyse called" << endl; - - if(fEventCounter==0) { + AliInfo("Analyse called"); - cout << "EventCounter == 0, Cannot analyse" << endl; + if(fIsAnalysed == kTRUE) { + + AliInfo("No new data since Analyse was called last time"); return; } + if(fEventCounter==0) { + + AliInfo("EventCounter == 0, Cannot analyse"); + return; + } + Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1; + AliInfo(Form("EventCounter: %d , TimeBins: %d", fEventCounter, nTimeBins)); + + Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins); + fNoThreshold->Multiply(normalization); - cout << "EventCounter: " << fEventCounter << endl - << "TimeBins: " << nTimeBins << endl; + fMeanCharge->Divide(fNLocalMaxima); + fMaxCharge->Divide(fNLocalMaxima); + fNTimeBins->Divide(fNLocalMaxima); + fNPads->Divide(fNLocalMaxima); + fTimePosition->Divide(fNLocalMaxima); + + fIsAnalysed = kTRUE; +} - if (fMeanCharge && fNoThreshold) fMeanCharge->Divide(fNoThreshold); - Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins); - if (fNoThreshold) fNoThreshold->Multiply(normalization); - if (fOverThreshold0) fOverThreshold0->Multiply(normalization); - if (fOverThreshold5) fOverThreshold5->Multiply(normalization); - if (fOverThreshold10) fOverThreshold10->Multiply(normalization); - if (fOverThreshold20) fOverThreshold20->Multiply(normalization); - if (fOverThreshold30) fOverThreshold30->Multiply(normalization); +//_____________________________________________________________________ +void AliTPCdataQA::MakeTree(const char *fname) const { + // + // Export result to the tree -located in the file + // This file can be analyzed using AliTPCCalibViewer + // + AliTPCPreprocessorOnline preprocesor; + + if (fNLocalMaxima) preprocesor.AddComponent(fNLocalMaxima); + if (fMaxCharge) preprocesor.AddComponent(fMaxCharge); + if (fMeanCharge) preprocesor.AddComponent(fMeanCharge); + if (fNoThreshold) preprocesor.AddComponent(fNoThreshold); + if (fNTimeBins) preprocesor.AddComponent(fNTimeBins); + if (fNPads) preprocesor.AddComponent(fNPads); + if (fTimePosition) preprocesor.AddComponent(fTimePosition); + if (fOverThreshold10) preprocesor.AddComponent(fOverThreshold10); + if (fOverThreshold20) preprocesor.AddComponent(fOverThreshold20); + if (fOverThreshold30) preprocesor.AddComponent(fOverThreshold30); + + preprocesor.DumpToFile(fname); +} + + +//_____________________________________________________________________ +void AliTPCdataQA::MakeArrays(){ + // + // The arrays for expanding the raw data are defined and + // som parameters are intialised + // + AliTPCROC * roc = AliTPCROC::Instance(); + // + // To make the array big enough for all sectors we take + // the dimensions from the outer row of an OROC (the last sector) + // + fRowsMax = roc->GetNRows(roc->GetNSector()-1); + fPadsMax = roc->GetNPads(roc->GetNSector()-1,fRowsMax-1); + fTimeBinsMax = fLastTimeBin - fFirstTimeBin +1; + + // + // Since we have added 2 pads (TimeBins) before and after the real pads (TimeBins) + // to make sure that we can always query the exanded table even when the + // max is on the edge + // + + + fAllBins = new Float_t*[fRowsMax]; + fAllSigBins = new Int_t*[fRowsMax]; + fAllNSigBins = new Int_t[fRowsMax]; + + for (Int_t iRow = 0; iRow < fRowsMax; iRow++) { + // + Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4); + fAllBins[iRow] = new Float_t[maxBin]; + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); // set all values to 0 + fAllSigBins[iRow] = new Int_t[maxBin]; + fAllNSigBins[iRow] = 0; + } +} + + +//_____________________________________________________________________ +void AliTPCdataQA::CleanArrays(){ + // + // + // + + for (Int_t iRow = 0; iRow < fRowsMax; iRow++) { + + // To speed up the performance by a factor 2 on cosmic data (and + // presumably pp data as well) where the ocupancy is low, the + // memset is only called if there is more than 1000 signals for a + // row (of the order 1% occupancy) + if(fAllNSigBins[iRow]<1000) { + + Float_t* allBins = fAllBins[iRow]; + Int_t* sigBins = fAllSigBins[iRow]; + const Int_t nSignals = fAllNSigBins[iRow]; + for(Int_t i = 0; i < nSignals; i++) + allBins[sigBins[i]]=0; + } else { + + Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4); + memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); + } + + fAllNSigBins[iRow]=0; + } +} + +//_____________________________________________________________________ +void AliTPCdataQA::GetPadAndTimeBin(Int_t bin, Int_t& iPad, Int_t& iTimeBin){ + // + // Return pad and timebin for a given bin + // + + // Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin; + iTimeBin = bin%(fTimeBinsMax+4); + iPad = (bin-iTimeBin)/(fTimeBinsMax+4); + + iPad -= 2; + iTimeBin -= 2; + iTimeBin += fFirstTimeBin; + + R__ASSERT(iPad>=0 && iPad<=fPadsMax); + R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin); +} + +//_____________________________________________________________________ +void AliTPCdataQA::SetExpandDigit(const Int_t iRow, Int_t iPad, + Int_t iTimeBin, const Float_t signal) +{ + // + // + // + R__ASSERT(iRow>=0 && iRow=0 && iPad<=fPadsMax); + R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin); + + iTimeBin -= fFirstTimeBin; + iPad += 2; + iTimeBin += 2; + + Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin; + fAllBins[iRow][bin] = signal; + fAllSigBins[iRow][fAllNSigBins[iRow]] = bin; + fAllNSigBins[iRow]++; +} + +Float_t AliTPCdataQA::GetQ(const Float_t* adcArray, const Int_t time, + const Int_t pad, const Int_t maxTimeBins, + Int_t& timeMin, Int_t& timeMax, + Int_t& padMin, Int_t& padMax) const +{ + // + // This methods return the charge in the bin time+pad*maxTimeBins + // If the charge is above 0 it also updates the padMin, padMax, timeMin + // and timeMax if necessary + // + Float_t charge = adcArray[time + pad*maxTimeBins]; + if(charge > 0) { + timeMin = TMath::Min(time, timeMin); timeMax = TMath::Max(time, timeMax); + padMin = TMath::Min(pad, padMin); padMax = TMath::Max(pad, padMax); + } + return charge; +} + +//______________________________________________________________________________ +void AliTPCdataQA::Streamer(TBuffer &xRuub) +{ + // Automatic schema evolution was first used from revision 4 + // Code based on: + // http://root.cern.ch/root/roottalk/roottalk02/3207.html + + UInt_t xRuus, xRuuc; + if (xRuub.IsReading()) { + Version_t xRuuv = xRuub.ReadVersion(&xRuus, &xRuuc); + //we use the automatic algorithm for class version > 3 + if (xRuuv > 3) { + AliTPCdataQA::Class()->ReadBuffer(xRuub, this, xRuuv, xRuus, + xRuuc); + return; + } + TH1F::Streamer(xRuub); + xRuub >> fFirstTimeBin; + xRuub >> fLastTimeBin; + xRuub >> fAdcMin; + xRuub >> fAdcMax; + xRuub >> fNLocalMaxima; + xRuub >> fMaxCharge; + xRuub >> fMeanCharge; + xRuub >> fNoThreshold; + xRuub >> fNTimeBins; + xRuub >> fNPads; + xRuub >> fTimePosition; + xRuub >> fEventCounter; + xRuub >> fIsAnalysed; + xRuub.CheckByteCount(xRuus, xRuuc, AliTPCdataQA::IsA()); + } else { + AliTPCdataQA::Class()->WriteBuffer(xRuub,this); + } }