1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
22 The code has been heavily modified so that now the RAW data is
23 "expanded" for each sector and stored in a big signal array. Then a
24 simple version of the code in AliTPCclustererMI is used to identify
25 the local maxima and these are then used for QA. This gives a better
26 estimate of the charge (both max and total) and also limits the
31 In Update the RAW signals >= 3 ADC channels are stored in the arrays.
34 Float_t** fAllBins 2d array [row][bin(pad, time)] ADC signal
35 Int_t** fAllSigBins 2d array [row][signal#] bin(with signal)
36 Int_t* fAllNSigBins; 1d array [row] Nsignals
38 This is done sector by sector.
40 When data from a new sector is encountered, the method
41 FindLocalMaxima is called on the data from the previous sector, and
42 the calibration/data objects are updated with the "cluster"
43 info. Finally the arrays are cleared.
45 The requirements for a local maxima is:
46 Charge in bin is >= 5 ADC channels.
47 Charge in bin is larger than all the 8 neighboring bins.
48 At least one of the two pad neighbors has a signal.
49 At least one of the two time neighbors has a signal.
51 Before accessing the data it is expected that the Analyse method is
52 called. This normalizes some of the data objects to per event or per
54 If more data is passed to the class after Analyse has been called
55 the normalization is reversed and Analyse has to be called again.
71 #include <TDirectory.h>
76 #include "AliRawReader.h"
77 #include "AliRawReaderRoot.h"
78 #include "AliRawReaderDate.h"
79 #include "AliTPCRawStream.h"
80 #include "AliTPCCalROC.h"
81 #include "AliTPCROC.h"
82 #include "AliMathBase.h"
83 #include "TTreeStream.h"
84 #include "AliTPCRawStreamFast.h"
88 #include "AliTPCCalPad.h"
89 #include "AliTPCPreprocessorOnline.h"
92 #include "AliTPCdataQA.h"
94 ClassImp(AliTPCdataQA)
96 AliTPCdataQA::AliTPCdataQA() : /*FOLD00*/
124 // default constructor
128 //_____________________________________________________________________
129 AliTPCdataQA::AliTPCdataQA(AliRecoParam::EventSpecie_t es) :
156 // ctor creating the histogram
157 char * name = Form("TPCRAW_%s", AliRecoParam::GetEventSpecieName(es)) ;
158 TH1F(name, name,100,0,100) ;
161 //_____________________________________________________________________
162 AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/
164 fFirstTimeBin(ped.GetFirstTimeBin()),
165 fLastTimeBin(ped.GetLastTimeBin()),
166 fAdcMin(ped.GetAdcMin()),
167 fAdcMax(ped.GetAdcMax()),
181 fEventCounter(ped.GetEventCounter()),
182 fIsAnalysed(ped.GetIsAnalysed()),
193 if(ped.GetNLocalMaxima())
194 fNLocalMaxima = new AliTPCCalPad(*ped.GetNLocalMaxima());
195 if(ped.GetMaxCharge())
196 fMaxCharge = new AliTPCCalPad(*ped.GetMaxCharge());
197 if(ped.GetMeanCharge())
198 fMeanCharge = new AliTPCCalPad(*ped.GetMeanCharge());
199 if(ped.GetNoThreshold())
200 fNoThreshold = new AliTPCCalPad(*ped.GetNoThreshold());
201 if(ped.GetNTimeBins())
202 fNTimeBins = new AliTPCCalPad(*ped.GetNTimeBins());
204 fNPads = new AliTPCCalPad(*ped.GetNPads());
205 if(ped.GetTimePosition())
206 fTimePosition = new AliTPCCalPad(*ped.GetTimePosition());
207 if(ped.GetOverThreshold10())
208 fOverThreshold10 = new AliTPCCalPad(*ped.GetOverThreshold10());
209 if(ped.GetOverThreshold20())
210 fOverThreshold20 = new AliTPCCalPad(*ped.GetOverThreshold20());
211 if(ped.GetOverThreshold30())
212 fOverThreshold30 = new AliTPCCalPad(*ped.GetOverThreshold30());
215 //_____________________________________________________________________
216 AliTPCdataQA::AliTPCdataQA(const TMap *config) : /*FOLD00*/
217 TH1F("TPCRAW","TPCRAW",100,0,100),
245 // default constructor
247 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
248 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
249 if (config->GetValue("AdcMin")) fAdcMin = ((TObjString*)config->GetValue("AdcMin"))->GetString().Atoi();
250 if (config->GetValue("AdcMax")) fAdcMax = ((TObjString*)config->GetValue("AdcMax"))->GetString().Atoi();
253 //_____________________________________________________________________
254 AliTPCdataQA& AliTPCdataQA::operator = (const AliTPCdataQA &source)
257 // assignment operator
259 if (&source == this) return *this;
260 new (this) AliTPCdataQA(source);
266 //_____________________________________________________________________
267 AliTPCdataQA::~AliTPCdataQA() /*FOLD00*/
273 // do not delete fMapping, because we do not own it.
274 // do not delete fMapping, because we do not own it.
275 // do not delete fNoise and fPedestal, because we do not own them.
277 delete fNLocalMaxima;
283 delete fTimePosition;
284 delete fOverThreshold10;
285 delete fOverThreshold20;
286 delete fOverThreshold30;
288 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
289 delete [] fAllBins[iRow];
290 delete [] fAllSigBins[iRow];
293 delete [] fAllSigBins;
294 delete [] fAllNSigBins;
297 //_____________________________________________________________________
298 Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
301 // Event Processing loop - AliTPCRawStream
303 Bool_t withInput = kFALSE;
305 Int_t lastSector = -1;
307 while ( rawStreamFast->NextDDL() ){
308 while ( rawStreamFast->NextChannel() ){
310 Int_t iSector = rawStreamFast->GetSector(); // current sector
311 Int_t iRow = rawStreamFast->GetRow(); // current row
312 Int_t iPad = rawStreamFast->GetPad(); // current pad
313 // Call local maxima finder if the data is in a new sector
314 if(iSector != lastSector) {
317 FindLocalMaxima(lastSector);
320 lastSector = iSector;
324 while ( rawStreamFast->NextBunch() ){
325 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
326 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
328 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
329 Float_t signal = rawStreamFast->GetSignals()[iTimeBin-startTbin];
330 nSignals += Update(iSector,iRow,iPad,iTimeBin+1,signal);
339 //_____________________________________________________________________
340 Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *rawReader)
343 // Event processing loop - AliRawReader
345 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
346 Bool_t res=ProcessEventFast(rawStreamFast);
348 fEventCounter++; // only increment event counter if there is TPC data
349 // otherwise Analyse (called in QA) fails
351 delete rawStreamFast;
355 //_____________________________________________________________________
356 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream)
359 // Event Processing loop - AliTPCRawStream
362 Bool_t withInput = kFALSE;
364 Int_t lastSector = -1;
366 while (rawStream->Next()) {
368 Int_t iSector = rawStream->GetSector(); // current ROC
369 Int_t iRow = rawStream->GetRow(); // current row
370 Int_t iPad = rawStream->GetPad(); // current pad
371 Int_t iTimeBin = rawStream->GetTime(); // current time bin
372 Float_t signal = rawStream->GetSignal(); // current ADC signal
374 // Call local maxima finder if the data is in a new sector
375 if(iSector != lastSector) {
378 FindLocalMaxima(lastSector);
381 lastSector = iSector;
385 // Sometimes iRow==-1 if there is problems to read the data
387 nSignals += Update(iSector,iRow,iPad,iTimeBin,signal);
396 //_____________________________________________________________________
397 Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader)
400 // Event processing loop - AliRawReader
403 // if fMapping is NULL the rawstream will crate its own mapping
404 AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
405 rawReader->Select("TPC");
406 Bool_t res = ProcessEvent(&rawStream);
409 fEventCounter++; // only increment event counter if there is TPC data
410 // otherwise Analyse (called in QA) fails
415 //_____________________________________________________________________
416 Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *event)
419 // process date event
422 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
423 Bool_t result=ProcessEvent(rawReader);
430 //_____________________________________________________________________
431 void AliTPCdataQA::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
434 // Write class to file
445 TDirectory *backup = gDirectory;
446 TFile f(filename,option.Data());
448 if ( !sDir.IsNull() ){
449 f.mkdir(sDir.Data());
455 if ( backup ) backup->cd();
459 //_____________________________________________________________________
460 Int_t AliTPCdataQA::Update(const Int_t iSector, /*FOLD00*/
463 const Int_t iTimeBin,
467 // Signal filling method
471 // Define the calibration objects the first time Update is called
472 // NB! This has to be done first even if the data is rejected by the time
473 // cut to make sure that the objects are available in Analyse
475 if (!fNLocalMaxima) fNLocalMaxima = new AliTPCCalPad("NLocalMaxima","NLocalMaxima");
476 if (!fMaxCharge) fMaxCharge = new AliTPCCalPad("MaxCharge","MaxCharge");
477 if (!fMeanCharge) fMeanCharge = new AliTPCCalPad("MeanCharge","MeanCharge");
478 if (!fNoThreshold) fNoThreshold = new AliTPCCalPad("NoThreshold","NoThreshold");
479 if (!fNTimeBins) fNTimeBins = new AliTPCCalPad("NTimeBins","NTimeBins");
480 if (!fNPads) fNPads = new AliTPCCalPad("NPads","NPads");
481 if (!fTimePosition) fTimePosition = new AliTPCCalPad("TimePosition","TimePosition");
482 if (!fOverThreshold10) fOverThreshold10 = new AliTPCCalPad("OverThreshold10","OverThreshold10");
483 if (!fOverThreshold20) fOverThreshold20 = new AliTPCCalPad("OverThreshold20","OverThreshold20");
484 if (!fOverThreshold30) fOverThreshold30 = new AliTPCCalPad("OverThreshold30","OverThreshold30");
485 // Make the arrays for expanding the data
491 // If Analyse has been previously called we need now to denormalize the data
492 // as more data is coming
494 if(fIsAnalysed == kTRUE) {
496 const Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
497 const Float_t denormalization = Float_t(fEventCounter * nTimeBins);
498 fNoThreshold->Multiply(denormalization);
500 fMeanCharge->Multiply(fNLocalMaxima);
501 fMaxCharge->Multiply(fNLocalMaxima);
502 fNTimeBins->Multiply(fNLocalMaxima);
503 fNPads->Multiply(fNLocalMaxima);
504 fTimePosition->Multiply(fNLocalMaxima);
505 fIsAnalysed = kFALSE;
511 if (iTimeBin<fFirstTimeBin) return 0;
512 if (iTimeBin>fLastTimeBin) return 0;
514 // if pedestal calibrations are loaded subtract pedestals
517 Float_t ped = fPedestal->GetCalROC(iSector)->GetValue(iRow, iPad);
518 // Don't use data from pads where pedestals are abnormally small or large
524 // In fNoThreshold we fill all data to estimate the ZS volume
525 Float_t count = fNoThreshold->GetCalROC(iSector)->GetValue(iRow, iPad);
526 fNoThreshold->GetCalROC(iSector)->SetValue(iRow, iPad,count+1);
528 // Require at least 3 ADC channels
532 // if noise calibrations are loaded require at least 3*sigmaNoise
535 Float_t noise = fNoise->GetCalROC(iSector)->GetValue(iRow, iPad);
537 if(signal < noise*3.0)
542 // This signal is ok and we store it in the cluster map
545 SetExpandDigit(iRow, iPad, iTimeBin, signal);
547 return 1; // signal was accepted
550 //_____________________________________________________________________
551 void AliTPCdataQA::FindLocalMaxima(const Int_t iSector)
554 // This method is called after the data from each sector has been
555 // exapanded into an array
556 // Loop over the signals and identify local maxima and fill the
557 // calibration objects with the information
560 Int_t nLocalMaxima = 0;
561 const Int_t maxTimeBin = fTimeBinsMax+4; // Used to step between neighboring pads
562 // Because we have tha pad-time data in a
565 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
567 Float_t* allBins = fAllBins[iRow];
568 Int_t* sigBins = fAllSigBins[iRow];
569 const Int_t nSigBins = fAllNSigBins[iRow];
571 for (Int_t iSig = 0; iSig < nSigBins; iSig++) {
573 Int_t bin = sigBins[iSig];
574 Float_t *b = &allBins[bin];
577 // Now we check if this is a local maximum
582 // First check that the charge is bigger than the threshold
586 // Require at least one neighboring pad with signal
587 if (b[-maxTimeBin]+b[maxTimeBin]<=0) continue;
589 // Require at least one neighboring time bin with signal
590 if (b[-1]+b[1]<=0) continue;
593 // Check that this is a local maximum
594 // Note that the checking is done so that if 2 charges has the same
595 // qMax then only 1 cluster is generated
596 // (that is why there is BOTH > and >=)
598 if (b[-maxTimeBin] >= qMax) continue;
599 if (b[-1 ] >= qMax) continue;
600 if (b[+maxTimeBin] > qMax) continue;
601 if (b[+1 ] > qMax) continue;
602 if (b[-maxTimeBin-1] >= qMax) continue;
603 if (b[+maxTimeBin-1] >= qMax) continue;
604 if (b[+maxTimeBin+1] > qMax) continue;
605 if (b[-maxTimeBin+1] >= qMax) continue;
608 // Now we accept the local maximum and fill the calibration/data objects
612 Int_t iPad, iTimeBin;
613 GetPadAndTimeBin(bin, iPad, iTimeBin);
615 Float_t count = fNLocalMaxima->GetCalROC(iSector)->GetValue(iRow, iPad);
616 fNLocalMaxima->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
618 count = fTimePosition->GetCalROC(iSector)->GetValue(iRow, iPad);
619 fTimePosition->GetCalROC(iSector)->SetValue(iRow, iPad, count+iTimeBin);
621 Float_t charge = fMaxCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
622 fMaxCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qMax);
625 count = fOverThreshold10->GetCalROC(iSector)->GetValue(iRow, iPad);
626 fOverThreshold10->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
629 count = fOverThreshold20->GetCalROC(iSector)->GetValue(iRow, iPad);
630 fOverThreshold20->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
633 count = fOverThreshold30->GetCalROC(iSector)->GetValue(iRow, iPad);
634 fOverThreshold30->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
638 // Calculate the total charge as the sum over the region:
646 // with qmax at the center C.
648 // The inner charge (i) we always add, but we only add the outer
649 // charge (o) if the neighboring inner bin (i) has a signal.
651 Int_t minP = 0, maxP = 0, minT = 0, maxT = 0;
653 for(Int_t i = -1; i<=1; i++) {
654 for(Int_t j = -1; j<=1; j++) {
659 Float_t charge1 = GetQ(b, i, j, maxTimeBin, minT, maxT, minP, maxP);
662 // see if the next neighbor is also above threshold
664 qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
666 // we are in a diagonal corner
667 qTot += GetQ(b, i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
668 qTot += GetQ(b, 2*i, j, maxTimeBin, minT, maxT, minP, maxP);
669 qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
675 charge = fMeanCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
676 fMeanCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qTot);
678 count = fNTimeBins->GetCalROC(iSector)->GetValue(iRow, iPad);
679 fNTimeBins->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxT-minT+1);
681 count = fNPads->GetCalROC(iSector)->GetValue(iRow, iPad);
682 fNPads->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxP-minP+1);
684 } // end loop over signals
685 } // end loop over rows
687 // cout << "Number of local maximas found: " << nLocalMaxima << endl;
690 //_____________________________________________________________________
691 void AliTPCdataQA::Analyse()
694 // Calculate calibration constants
697 cout << "Analyse called" << endl;
699 if(fIsAnalysed == kTRUE) {
701 cout << "No new data since Analyse was called last time" << endl;
705 if(fEventCounter==0) {
707 cout << "EventCounter == 0, Cannot analyse" << endl;
711 Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
712 cout << "EventCounter: " << fEventCounter << endl
713 << "TimeBins: " << nTimeBins << endl;
715 Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins);
716 fNoThreshold->Multiply(normalization);
718 fMeanCharge->Divide(fNLocalMaxima);
719 fMaxCharge->Divide(fNLocalMaxima);
720 fNTimeBins->Divide(fNLocalMaxima);
721 fNPads->Divide(fNLocalMaxima);
722 fTimePosition->Divide(fNLocalMaxima);
728 //_____________________________________________________________________
729 void AliTPCdataQA::MakeTree(const char *fname){
731 // Export result to the tree -located in the file
732 // This file can be analyzed using AliTPCCalibViewer
734 AliTPCPreprocessorOnline preprocesor;
736 if (fNLocalMaxima) preprocesor.AddComponent(fNLocalMaxima);
737 if (fMaxCharge) preprocesor.AddComponent(fMaxCharge);
738 if (fMeanCharge) preprocesor.AddComponent(fMeanCharge);
739 if (fNoThreshold) preprocesor.AddComponent(fNoThreshold);
740 if (fNTimeBins) preprocesor.AddComponent(fNTimeBins);
741 if (fNPads) preprocesor.AddComponent(fNPads);
742 if (fTimePosition) preprocesor.AddComponent(fTimePosition);
743 if (fOverThreshold10) preprocesor.AddComponent(fOverThreshold10);
744 if (fOverThreshold20) preprocesor.AddComponent(fOverThreshold20);
745 if (fOverThreshold30) preprocesor.AddComponent(fOverThreshold30);
747 preprocesor.DumpToFile(fname);
751 //_____________________________________________________________________
752 void AliTPCdataQA::MakeArrays(){
754 // The arrays for expanding the raw data are defined and
755 // som parameters are intialised
757 AliTPCROC * roc = AliTPCROC::Instance();
759 // To make the array big enough for all sectors we take
760 // the dimensions from the outer row of an OROC (the last sector)
762 fRowsMax = roc->GetNRows(roc->GetNSector()-1);
763 fPadsMax = roc->GetNPads(roc->GetNSector()-1,fRowsMax-1);
764 fTimeBinsMax = fLastTimeBin - fFirstTimeBin +1;
767 // Since we have added 2 pads (TimeBins) before and after the real pads (TimeBins)
768 // to make sure that we can always query the exanded table even when the
769 // max is on the edge
773 fAllBins = new Float_t*[fRowsMax];
774 fAllSigBins = new Int_t*[fRowsMax];
775 fAllNSigBins = new Int_t[fRowsMax];
777 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
779 Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4);
780 fAllBins[iRow] = new Float_t[maxBin];
781 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); // set all values to 0
782 fAllSigBins[iRow] = new Int_t[maxBin];
783 fAllNSigBins[iRow] = 0;
788 //_____________________________________________________________________
789 void AliTPCdataQA::CleanArrays(){
794 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
796 Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4);
797 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
798 fAllNSigBins[iRow]=0;
802 //_____________________________________________________________________
803 void AliTPCdataQA::GetPadAndTimeBin(Int_t bin, Int_t& iPad, Int_t& iTimeBin){
805 // Return pad and timebin for a given bin
808 // Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
809 iTimeBin = bin%(fTimeBinsMax+4);
810 iPad = (bin-iTimeBin)/(fTimeBinsMax+4);
814 iTimeBin += fFirstTimeBin;
816 R__ASSERT(iPad>=0 && iPad<=fPadsMax);
817 R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
820 //_____________________________________________________________________
821 void AliTPCdataQA::SetExpandDigit(const Int_t iRow, Int_t iPad,
822 Int_t iTimeBin, const Float_t signal)
827 R__ASSERT(iRow>=0 && iRow<fRowsMax);
828 R__ASSERT(iPad>=0 && iPad<=fPadsMax);
829 R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
831 iTimeBin -= fFirstTimeBin;
835 Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
836 fAllBins[iRow][bin] = signal;
837 fAllSigBins[iRow][fAllNSigBins[iRow]] = bin;
838 fAllNSigBins[iRow]++;
841 Float_t AliTPCdataQA::GetQ(const Float_t* adcArray, const Int_t time,
842 const Int_t pad, const Int_t maxTimeBins,
843 Int_t& timeMin, Int_t& timeMax,
844 Int_t& padMin, Int_t& padMax)
847 // This methods return the charge in the bin time+pad*maxTimeBins
848 // If the charge is above 0 it also updates the padMin, padMax, timeMin
849 // and timeMax if necessary
851 Float_t charge = adcArray[time + pad*maxTimeBins];
853 timeMin = TMath::Min(time, timeMin); timeMax = TMath::Max(time, timeMax);
854 padMin = TMath::Min(pad, padMin); padMax = TMath::Max(pad, padMax);