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>
75 #include "AliRawReader.h"
76 #include "AliRawReaderRoot.h"
77 #include "AliRawReaderDate.h"
78 #include "AliTPCRawStream.h"
79 #include "AliTPCCalROC.h"
80 #include "AliTPCROC.h"
81 #include "AliMathBase.h"
82 #include "TTreeStream.h"
83 #include "AliTPCRawStreamFast.h"
87 #include "AliTPCCalPad.h"
88 #include "AliTPCPreprocessorOnline.h"
91 #include "AliTPCdataQA.h"
93 ClassImp(AliTPCdataQA)
95 AliTPCdataQA::AliTPCdataQA() : /*FOLD00*/
96 TH1F("TPCRAW","TPCRAW",100,0,100),
101 fOldRCUformat(kTRUE),
125 // default constructor
130 //_____________________________________________________________________
131 AliTPCdataQA::AliTPCdataQA(const AliTPCdataQA &ped) : /*FOLD00*/
133 fFirstTimeBin(ped.GetFirstTimeBin()),
134 fLastTimeBin(ped.GetLastTimeBin()),
135 fAdcMin(ped.GetAdcMin()),
136 fAdcMax(ped.GetAdcMax()),
137 fOldRCUformat(ped.GetOldRCUformat()),
151 fEventCounter(ped.GetEventCounter()),
152 fIsAnalysed(ped.GetIsAnalysed()),
163 if(ped.GetNLocalMaxima())
164 fNLocalMaxima = new AliTPCCalPad(*ped.GetNLocalMaxima());
165 if(ped.GetMaxCharge())
166 fMaxCharge = new AliTPCCalPad(*ped.GetMaxCharge());
167 if(ped.GetMeanCharge())
168 fMeanCharge = new AliTPCCalPad(*ped.GetMeanCharge());
169 if(ped.GetNoThreshold())
170 fNoThreshold = new AliTPCCalPad(*ped.GetNoThreshold());
171 if(ped.GetNTimeBins())
172 fNTimeBins = new AliTPCCalPad(*ped.GetNTimeBins());
174 fNPads = new AliTPCCalPad(*ped.GetNPads());
175 if(ped.GetTimePosition())
176 fTimePosition = new AliTPCCalPad(*ped.GetTimePosition());
177 if(ped.GetOverThreshold10())
178 fOverThreshold10 = new AliTPCCalPad(*ped.GetOverThreshold10());
179 if(ped.GetOverThreshold20())
180 fOverThreshold20 = new AliTPCCalPad(*ped.GetOverThreshold20());
181 if(ped.GetOverThreshold30())
182 fOverThreshold30 = new AliTPCCalPad(*ped.GetOverThreshold30());
186 //_____________________________________________________________________
187 AliTPCdataQA& AliTPCdataQA::operator = (const AliTPCdataQA &source)
190 // assignment operator
192 if (&source == this) return *this;
193 new (this) AliTPCdataQA(source);
199 //_____________________________________________________________________
200 AliTPCdataQA::~AliTPCdataQA() /*FOLD00*/
206 // do not delete fMapping, because we do not own it.
207 // do not delete fMapping, because we do not own it.
208 // do not delete fNoise and fPedestal, because we do not own them.
210 delete fNLocalMaxima;
216 delete fTimePosition;
217 delete fOverThreshold10;
218 delete fOverThreshold20;
219 delete fOverThreshold30;
221 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
222 delete [] fAllBins[iRow];
223 delete [] fAllSigBins[iRow];
226 delete [] fAllSigBins;
227 delete [] fAllNSigBins;
230 //_____________________________________________________________________
231 Bool_t AliTPCdataQA::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
234 // Event Processing loop - AliTPCRawStream
236 Bool_t withInput = kFALSE;
238 Int_t lastSector = -1;
240 while ( rawStreamFast->NextDDL() ){
241 while ( rawStreamFast->NextChannel() ){
243 Int_t iSector = rawStreamFast->GetSector(); // current sector
244 Int_t iRow = rawStreamFast->GetRow(); // current row
245 Int_t iPad = rawStreamFast->GetPad(); // current pad
246 // Call local maxima finder if the data is in a new sector
247 if(iSector != lastSector) {
250 FindLocalMaxima(lastSector);
253 lastSector = iSector;
257 while ( rawStreamFast->NextBunch() ){
258 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
259 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
261 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
262 Float_t signal = rawStreamFast->GetSignals()[iTimeBin-startTbin];
263 nSignals += Update(iSector,iRow,iPad,iTimeBin+1,signal);
272 //_____________________________________________________________________
273 Bool_t AliTPCdataQA::ProcessEventFast(AliRawReader *rawReader)
276 // Event processing loop - AliRawReader
278 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
279 Bool_t res=ProcessEventFast(rawStreamFast);
281 fEventCounter++; // only increment event counter if there is TPC data
282 // otherwise Analyse (called in QA) fails
284 delete rawStreamFast;
288 //_____________________________________________________________________
289 Bool_t AliTPCdataQA::ProcessEvent(AliTPCRawStream *rawStream)
292 // Event Processing loop - AliTPCRawStream
295 rawStream->SetOldRCUFormat(fOldRCUformat);
297 Bool_t withInput = kFALSE;
299 Int_t lastSector = -1;
301 while (rawStream->Next()) {
303 Int_t iSector = rawStream->GetSector(); // current ROC
304 Int_t iRow = rawStream->GetRow(); // current row
305 Int_t iPad = rawStream->GetPad(); // current pad
306 Int_t iTimeBin = rawStream->GetTime(); // current time bin
307 Float_t signal = rawStream->GetSignal(); // current ADC signal
309 // Call local maxima finder if the data is in a new sector
310 if(iSector != lastSector) {
313 FindLocalMaxima(lastSector);
316 lastSector = iSector;
320 // Sometimes iRow==-1 if there is problems to read the data
322 nSignals += Update(iSector,iRow,iPad,iTimeBin,signal);
331 //_____________________________________________________________________
332 Bool_t AliTPCdataQA::ProcessEvent(AliRawReader *rawReader)
335 // Event processing loop - AliRawReader
338 // if fMapping is NULL the rawstream will crate its own mapping
339 AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
340 rawReader->Select("TPC");
341 Bool_t res = ProcessEvent(&rawStream);
344 fEventCounter++; // only increment event counter if there is TPC data
345 // otherwise Analyse (called in QA) fails
350 //_____________________________________________________________________
351 Bool_t AliTPCdataQA::ProcessEvent(eventHeaderStruct *event)
354 // process date event
357 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
358 Bool_t result=ProcessEvent(rawReader);
365 //_____________________________________________________________________
366 void AliTPCdataQA::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
369 // Write class to file
380 TDirectory *backup = gDirectory;
381 TFile f(filename,option.Data());
383 if ( !sDir.IsNull() ){
384 f.mkdir(sDir.Data());
390 if ( backup ) backup->cd();
394 //_____________________________________________________________________
395 Int_t AliTPCdataQA::Update(const Int_t iSector, /*FOLD00*/
398 const Int_t iTimeBin,
402 // Signal filling method
406 // Define the calibration objects the first time Update is called
407 // NB! This has to be done first even if the data is rejected by the time
408 // cut to make sure that the objects are available in Analyse
410 if (!fNLocalMaxima) fNLocalMaxima = new AliTPCCalPad("NLocalMaxima","NLocalMaxima");
411 if (!fMaxCharge) fMaxCharge = new AliTPCCalPad("MaxCharge","MaxCharge");
412 if (!fMeanCharge) fMeanCharge = new AliTPCCalPad("MeanCharge","MeanCharge");
413 if (!fNoThreshold) fNoThreshold = new AliTPCCalPad("NoThreshold","NoThreshold");
414 if (!fNTimeBins) fNTimeBins = new AliTPCCalPad("NTimeBins","NTimeBins");
415 if (!fNPads) fNPads = new AliTPCCalPad("NPads","NPads");
416 if (!fTimePosition) fTimePosition = new AliTPCCalPad("TimePosition","TimePosition");
417 if (!fOverThreshold10) fOverThreshold10 = new AliTPCCalPad("OverThreshold10","OverThreshold10");
418 if (!fOverThreshold20) fOverThreshold20 = new AliTPCCalPad("OverThreshold20","OverThreshold20");
419 if (!fOverThreshold30) fOverThreshold30 = new AliTPCCalPad("OverThreshold30","OverThreshold30");
420 // Make the arrays for expanding the data
426 // If Analyse has been previously called we need now to denormalize the data
427 // as more data is coming
429 if(fIsAnalysed == kTRUE) {
431 const Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
432 const Float_t denormalization = Float_t(fEventCounter * nTimeBins);
433 fNoThreshold->Multiply(denormalization);
435 fMeanCharge->Multiply(fNLocalMaxima);
436 fMaxCharge->Multiply(fNLocalMaxima);
437 fNTimeBins->Multiply(fNLocalMaxima);
438 fNPads->Multiply(fNLocalMaxima);
439 fTimePosition->Multiply(fNLocalMaxima);
440 fIsAnalysed = kFALSE;
446 if (iTimeBin<fFirstTimeBin) return 0;
447 if (iTimeBin>fLastTimeBin) return 0;
449 // if pedestal calibrations are loaded subtract pedestals
452 Float_t ped = fPedestal->GetCalROC(iSector)->GetValue(iRow, iPad);
453 // Don't use data from pads where pedestals are abnormally small or large
459 // In fNoThreshold we fill all data to estimate the ZS volume
460 Float_t count = fNoThreshold->GetCalROC(iSector)->GetValue(iRow, iPad);
461 fNoThreshold->GetCalROC(iSector)->SetValue(iRow, iPad,count+1);
463 // Require at least 3 ADC channels
467 // if noise calibrations are loaded require at least 3*sigmaNoise
470 Float_t noise = fNoise->GetCalROC(iSector)->GetValue(iRow, iPad);
472 if(signal < noise*3.0)
477 // This signal is ok and we store it in the cluster map
480 SetExpandDigit(iRow, iPad, iTimeBin, signal);
482 return 1; // signal was accepted
485 //_____________________________________________________________________
486 void AliTPCdataQA::FindLocalMaxima(const Int_t iSector)
489 // This method is called after the data from each sector has been
490 // exapanded into an array
491 // Loop over the signals and identify local maxima and fill the
492 // calibration objects with the information
495 Int_t nLocalMaxima = 0;
496 const Int_t maxTimeBin = fTimeBinsMax+4; // Used to step between neighboring pads
497 // Because we have tha pad-time data in a
500 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
502 Float_t* allBins = fAllBins[iRow];
503 Int_t* sigBins = fAllSigBins[iRow];
504 const Int_t nSigBins = fAllNSigBins[iRow];
506 for (Int_t iSig = 0; iSig < nSigBins; iSig++) {
508 Int_t bin = sigBins[iSig];
509 Float_t *b = &allBins[bin];
512 // Now we check if this is a local maximum
517 // First check that the charge is bigger than the threshold
521 // Require at least one neighboring pad with signal
522 if (b[-maxTimeBin]+b[maxTimeBin]<=0) continue;
524 // Require at least one neighboring time bin with signal
525 if (b[-1]+b[1]<=0) continue;
528 // Check that this is a local maximum
529 // Note that the checking is done so that if 2 charges has the same
530 // qMax then only 1 cluster is generated
531 // (that is why there is BOTH > and >=)
533 if (b[-maxTimeBin] >= qMax) continue;
534 if (b[-1 ] >= qMax) continue;
535 if (b[+maxTimeBin] > qMax) continue;
536 if (b[+1 ] > qMax) continue;
537 if (b[-maxTimeBin-1] >= qMax) continue;
538 if (b[+maxTimeBin-1] >= qMax) continue;
539 if (b[+maxTimeBin+1] > qMax) continue;
540 if (b[-maxTimeBin+1] >= qMax) continue;
543 // Now we accept the local maximum and fill the calibration/data objects
547 Int_t iPad, iTimeBin;
548 GetPadAndTimeBin(bin, iPad, iTimeBin);
550 Float_t count = fNLocalMaxima->GetCalROC(iSector)->GetValue(iRow, iPad);
551 fNLocalMaxima->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
553 count = fTimePosition->GetCalROC(iSector)->GetValue(iRow, iPad);
554 fTimePosition->GetCalROC(iSector)->SetValue(iRow, iPad, count+iTimeBin);
556 Float_t charge = fMaxCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
557 fMaxCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qMax);
560 count = fOverThreshold10->GetCalROC(iSector)->GetValue(iRow, iPad);
561 fOverThreshold10->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
564 count = fOverThreshold20->GetCalROC(iSector)->GetValue(iRow, iPad);
565 fOverThreshold20->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
568 count = fOverThreshold30->GetCalROC(iSector)->GetValue(iRow, iPad);
569 fOverThreshold30->GetCalROC(iSector)->SetValue(iRow, iPad, count+1);
573 // Calculate the total charge as the sum over the region:
581 // with qmax at the center C.
583 // The inner charge (i) we always add, but we only add the outer
584 // charge (o) if the neighboring inner bin (i) has a signal.
586 Int_t minP = 0, maxP = 0, minT = 0, maxT = 0;
588 for(Int_t i = -1; i<=1; i++) {
589 for(Int_t j = -1; j<=1; j++) {
594 Float_t charge = GetQ(b, i, j, maxTimeBin, minT, maxT, minP, maxP);
597 // see if the next neighbor is also above threshold
599 qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
601 // we are in a diagonal corner
602 qTot += GetQ(b, i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
603 qTot += GetQ(b, 2*i, j, maxTimeBin, minT, maxT, minP, maxP);
604 qTot += GetQ(b, 2*i, 2*j, maxTimeBin, minT, maxT, minP, maxP);
610 charge = fMeanCharge->GetCalROC(iSector)->GetValue(iRow, iPad);
611 fMeanCharge->GetCalROC(iSector)->SetValue(iRow, iPad, charge + qTot);
613 count = fNTimeBins->GetCalROC(iSector)->GetValue(iRow, iPad);
614 fNTimeBins->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxT-minT+1);
616 count = fNPads->GetCalROC(iSector)->GetValue(iRow, iPad);
617 fNPads->GetCalROC(iSector)->SetValue(iRow, iPad, count + maxP-minP+1);
619 } // end loop over signals
620 } // end loop over rows
622 // cout << "Number of local maximas found: " << nLocalMaxima << endl;
625 //_____________________________________________________________________
626 void AliTPCdataQA::Analyse()
629 // Calculate calibration constants
632 cout << "Analyse called" << endl;
634 if(fIsAnalysed == kTRUE) {
636 cout << "No new data since Analyse was called last time" << endl;
640 if(fEventCounter==0) {
642 cout << "EventCounter == 0, Cannot analyse" << endl;
646 Int_t nTimeBins = fLastTimeBin - fFirstTimeBin +1;
647 cout << "EventCounter: " << fEventCounter << endl
648 << "TimeBins: " << nTimeBins << endl;
650 Float_t normalization = 1.0 / Float_t(fEventCounter * nTimeBins);
651 fNoThreshold->Multiply(normalization);
653 fMeanCharge->Divide(fNLocalMaxima);
654 fMaxCharge->Divide(fNLocalMaxima);
655 fNTimeBins->Divide(fNLocalMaxima);
656 fNPads->Divide(fNLocalMaxima);
657 fTimePosition->Divide(fNLocalMaxima);
663 //_____________________________________________________________________
664 void AliTPCdataQA::MakeTree(const char *fname){
666 // Export result to the tree -located in the file
667 // This file can be analyzed using AliTPCCalibViewer
669 AliTPCPreprocessorOnline preprocesor;
671 if (fNLocalMaxima) preprocesor.AddComponent(fNLocalMaxima);
672 if (fMaxCharge) preprocesor.AddComponent(fMaxCharge);
673 if (fMeanCharge) preprocesor.AddComponent(fMeanCharge);
674 if (fNoThreshold) preprocesor.AddComponent(fNoThreshold);
675 if (fNTimeBins) preprocesor.AddComponent(fNTimeBins);
676 if (fNPads) preprocesor.AddComponent(fNPads);
677 if (fTimePosition) preprocesor.AddComponent(fTimePosition);
678 if (fOverThreshold10) preprocesor.AddComponent(fOverThreshold10);
679 if (fOverThreshold20) preprocesor.AddComponent(fOverThreshold20);
680 if (fOverThreshold30) preprocesor.AddComponent(fOverThreshold30);
682 preprocesor.DumpToFile(fname);
686 //_____________________________________________________________________
687 void AliTPCdataQA::MakeArrays(){
689 // The arrays for expanding the raw data are defined and
690 // som parameters are intialised
692 AliTPCROC * roc = AliTPCROC::Instance();
694 // To make the array big enough for all sectors we take
695 // the dimensions from the outer row of an OROC (the last sector)
697 fRowsMax = roc->GetNRows(roc->GetNSector()-1);
698 fPadsMax = roc->GetNPads(roc->GetNSector()-1,fRowsMax-1);
699 fTimeBinsMax = fLastTimeBin - fFirstTimeBin +1;
702 // Since we have added 2 pads (TimeBins) before and after the real pads (TimeBins)
703 // to make sure that we can always query the exanded table even when the
704 // max is on the edge
708 fAllBins = new Float_t*[fRowsMax];
709 fAllSigBins = new Int_t*[fRowsMax];
710 fAllNSigBins = new Int_t[fRowsMax];
712 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
714 Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4);
715 fAllBins[iRow] = new Float_t[maxBin];
716 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin); // set all values to 0
717 fAllSigBins[iRow] = new Int_t[maxBin];
718 fAllNSigBins[iRow] = 0;
723 //_____________________________________________________________________
724 void AliTPCdataQA::CleanArrays(){
729 for (Int_t iRow = 0; iRow < fRowsMax; iRow++) {
731 Int_t maxBin = (fTimeBinsMax+4)*(fPadsMax+4);
732 memset(fAllBins[iRow],0,sizeof(Float_t)*maxBin);
733 fAllNSigBins[iRow]=0;
737 //_____________________________________________________________________
738 void AliTPCdataQA::GetPadAndTimeBin(Int_t bin, Int_t& iPad, Int_t& iTimeBin){
740 // Return pad and timebin for a given bin
743 // Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
744 iTimeBin = bin%(fTimeBinsMax+4);
745 iPad = (bin-iTimeBin)/(fTimeBinsMax+4);
749 iTimeBin += fFirstTimeBin;
751 R__ASSERT(iPad>=0 && iPad<=fPadsMax);
752 R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
755 //_____________________________________________________________________
756 void AliTPCdataQA::SetExpandDigit(const Int_t iRow, Int_t iPad,
757 Int_t iTimeBin, const Float_t signal)
762 R__ASSERT(iRow>=0 && iRow<fRowsMax);
763 R__ASSERT(iPad>=0 && iPad<=fPadsMax);
764 R__ASSERT(iTimeBin>=fFirstTimeBin && iTimeBin<=fLastTimeBin);
766 iTimeBin -= fFirstTimeBin;
770 Int_t bin = iPad*(fTimeBinsMax+4)+iTimeBin;
771 fAllBins[iRow][bin] = signal;
772 fAllSigBins[iRow][fAllNSigBins[iRow]] = bin;
773 fAllNSigBins[iRow]++;
776 Float_t AliTPCdataQA::GetQ(const Float_t* adcArray, const Int_t time,
777 const Int_t pad, const Int_t maxTimeBins,
778 Int_t& timeMin, Int_t& timeMax,
779 Int_t& padMin, Int_t& padMax)
782 // This methods return the charge in the bin time+pad*maxTimeBins
783 // If the charge is above 0 it also updates the padMin, padMax, timeMin
784 // and timeMax if necessary
786 Float_t charge = adcArray[time + pad*maxTimeBins];
788 timeMin = TMath::Min(time, timeMin); timeMax = TMath::Max(time, timeMax);
789 padMin = TMath::Min(pad, padMin); padMax = TMath::Max(pad, padMax);