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 **************************************************************************/
19 example: fill pedestal with Gaussian noise
20 AliTRDCalibPadStatus ped;
21 ped.TestEvent(numberofevent);
22 // Method without histo
23 //////////////////////////
25 //Create the histo of the AliTRDCalROC
26 TH2F * histo2dm = ped.GetCalRocMean(0,kFALSE)->MakeHisto2D();
27 histo2dm->Scale(10.0);
28 TH1F * histo1dm = ped.GetCalRocMean(0,kFALSE)->MakeHisto1D();
29 histo1dm->Scale(10.0);
30 TH2F * histo2ds = ped.GetCalRocSquares(0,kFALSE)->MakeHisto2D();
31 histo2ds->Scale(10.0);
32 TH1F * histo1ds = ped.GetCalRocSquares(0,kFALSE)->MakeHisto1D();
35 TCanvas* c1 = new TCanvas;
38 histo2dm->Draw("colz");
42 histo2ds->Draw("colz");
46 /////////////////////////
49 TH1F * histo = ped.GetHisto(31);
56 #include <TObjArray.h>
63 #include <TDirectory.h>
65 #include <TTreeStream.h>
68 #include <AliMathBase.h>
69 #include "AliRawReader.h"
70 #include "AliRawReaderRoot.h"
71 #include "AliRawReaderDate.h"
74 #include "AliTRDCalibPadStatus.h"
75 #include "AliTRDRawStream.h"
76 #include "AliTRDarrayF.h"
77 #include "AliTRDgeometry.h"
78 #include "AliTRDCommonParam.h"
79 #include "./Cal/AliTRDCalROC.h"
80 #include "./Cal/AliTRDCalPadStatus.h"
81 #include "./Cal/AliTRDCalSingleChamberStatus.h"
87 ClassImp(AliTRDCalibPadStatus) /*FOLD00*/
89 //_____________________________________________________________________
90 AliTRDCalibPadStatus::AliTRDCalibPadStatus() : /*FOLD00*/
97 fCalArrayEntries(540),
99 fCalArraySquares(540),
100 fCalRocArrayMean(540),
101 fCalRocArrayRMS(540),
108 // default constructor
111 fGeo = new AliTRDgeometry();
115 //_____________________________________________________________________
116 AliTRDCalibPadStatus::AliTRDCalibPadStatus(const AliTRDCalibPadStatus &ped) : /*FOLD00*/
119 fAdcMin(ped.GetAdcMin()),
120 fAdcMax(ped.GetAdcMax()),
121 fDetector(ped.fDetector),
122 fNumberOfTimeBins(ped.fNumberOfTimeBins),
123 fCalArrayEntries(540),
125 fCalArraySquares(540),
126 fCalRocArrayMean(540),
127 fCalRocArrayRMS(540),
136 for (Int_t idet = 0; idet < 540; idet++){
137 const AliTRDarrayF *calEntries = (AliTRDarrayF*)ped.fCalArrayEntries.UncheckedAt(idet);
138 const AliTRDarrayF *calMean = (AliTRDarrayF*)ped.fCalArrayMean.UncheckedAt(idet);
139 const AliTRDarrayF *calSquares = (AliTRDarrayF*)ped.fCalArraySquares.UncheckedAt(idet);
140 const AliTRDCalROC *calRocMean = (AliTRDCalROC*)ped.fCalRocArrayMean.UncheckedAt(idet);
141 const AliTRDCalROC *calRocRMS = (AliTRDCalROC*)ped.fCalRocArrayRMS.UncheckedAt(idet);
142 const TH2F *hped = (TH2F*)ped.fHistoArray.UncheckedAt(idet);
144 if ( calEntries != 0x0 ) fCalArrayEntries.AddAt(new AliTRDarrayF(*calEntries), idet);
145 if ( calMean != 0x0 ) fCalArrayMean.AddAt(new AliTRDarrayF(*calMean), idet);
146 if ( calSquares != 0x0 ) fCalArraySquares.AddAt(new AliTRDarrayF(*calSquares), idet);
147 if ( calRocMean != 0x0 ) fCalRocArrayMean.AddAt(new AliTRDCalROC(*calRocMean), idet);
148 if ( calRocRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTRDCalROC(*calRocRMS), idet);
151 TH2F *hNew = new TH2F(*hped);
152 hNew->SetDirectory(0);
153 fHistoArray.AddAt(hNew,idet);
160 fGeo = new AliTRDgeometry();
163 //_____________________________________________________________________
164 AliTRDCalibPadStatus& AliTRDCalibPadStatus::operator = (const AliTRDCalibPadStatus &source)
167 // assignment operator
169 if (&source == this) return *this;
170 new (this) AliTRDCalibPadStatus(source);
175 //_____________________________________________________________________
176 AliTRDCalibPadStatus::~AliTRDCalibPadStatus() /*FOLD00*/
186 //_____________________________________________________________________
187 Int_t AliTRDCalibPadStatus::Update(const Int_t icdet, /*FOLD00*/
194 // Signal filling methode
196 if ( (csignal>fAdcMax) || (csignal<fAdcMin) ) return 0;
198 if(fDetector != icdet){
199 fCalEntries = ((AliTRDarrayF *)GetCalEntries(icdet,kTRUE));
200 fCalMean = ((AliTRDarrayF *)GetCalMean(icdet,kTRUE));
201 fCalSquares = ((AliTRDarrayF *)GetCalSquares(icdet,kTRUE));
204 Float_t entries = fCalEntries->At(icRow+icCol*crowMax);
205 Float_t mean = fCalMean->At(icRow+icCol*crowMax);
206 Float_t squares = fCalSquares->At(icRow+icCol*crowMax);
208 Float_t entriesn = entries+1.0;
209 fCalEntries->AddAt(entriesn,(icRow+icCol*crowMax));
210 Float_t meann = (mean*entries+((Float_t)(csignal+0.5)))/entriesn;
211 fCalMean->AddAt(meann,icRow+icCol*crowMax);
212 Float_t squaresn = ((squares*entries)+(((Float_t)(csignal+0.5))*((Float_t)(csignal+0.5))))/entriesn;
213 fCalSquares->AddAt(squaresn,icRow+icCol*crowMax);
215 //printf("icdet %d, icRow %d, icCol %d, csignal %d, crowMax %d\n",icdet,icRow,icCol,csignal,crowMax);
216 //printf("entries %f, mean %f, squares %f\n",entriesn,meann,squaresn);
223 //_____________________________________________________________________
224 Int_t AliTRDCalibPadStatus::UpdateHisto(const Int_t icdet, /*FOLD00*/
231 // Signal filling methode
233 Int_t nbchannel = icRow+icCol*crowMax;
235 // fast filling methode.
236 // Attention: the entry counter of the histogram is not increased
237 // this means that e.g. the colz draw option gives an empty plot
239 if ( !(((Int_t)csignal>fAdcMax ) || ((Int_t)csignal<fAdcMin)) )
240 bin = (nbchannel+1)*(fAdcMax-fAdcMin+2)+((Int_t)csignal-fAdcMin+1);
242 GetHisto(icdet,kTRUE)->GetArray()[bin]++;
247 //_____________________________________________________________________
248 Int_t AliTRDCalibPadStatus::ProcessEvent(AliTRDRawStream *rawStream, Bool_t nocheck)
251 // Event Processing loop - AliTRDRawStream
252 // 0 time bin problem or zero suppression
260 while (rawStream->Next()) {
261 Int_t rawversion = rawStream->GetRawVersion(); // current raw version
262 if(rawversion != 2) return 0;
263 Int_t idetector = rawStream->GetDet(); // current detector
264 Int_t iRow = rawStream->GetRow(); // current row
265 Int_t iRowMax = rawStream->GetMaxRow(); // current rowmax
266 Int_t iCol = rawStream->GetCol(); // current col
267 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
268 Int_t *signal = rawStream->GetSignals(); // current ADC signal
269 Int_t nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
271 if((fDetector != -1) && (nbtimebin != fNumberOfTimeBins)) return 0;
272 fNumberOfTimeBins = nbtimebin;
274 Int_t fin = TMath::Min(nbtimebin,(iTimeBin+3));
277 for(Int_t k = iTimeBin; k < fin; k++){
278 if(signal[n]>0) UpdateHisto(idetector,iRow,iCol,signal[n],iRowMax);
286 while (rawStream->Next()) {
287 Int_t idetector = rawStream->GetDet(); // current detector
288 Int_t iRow = rawStream->GetRow(); // current row
289 Int_t iRowMax = rawStream->GetMaxRow(); // current rowmax
290 Int_t iCol = rawStream->GetCol(); // current col
291 Int_t iTimeBin = rawStream->GetTimeBin(); // current time bin
292 Int_t *signal = rawStream->GetSignals(); // current ADC signal
293 Int_t nbtimebin = rawStream->GetNumberOfTimeBins(); // number of time bins read from data
295 Int_t fin = TMath::Min(nbtimebin,(iTimeBin+3));
298 for(Int_t k = iTimeBin; k < fin; k++){
299 if(signal[n]>0) UpdateHisto(idetector,iRow,iCol,signal[n],iRowMax);
310 //_____________________________________________________________________
311 Int_t AliTRDCalibPadStatus::ProcessEvent(AliRawReader *rawReader, Bool_t nocheck)
314 // Event processing loop - AliRawReader
318 AliTRDRawStream rawStream(rawReader);
320 rawReader->Select("TRD");
322 return ProcessEvent(&rawStream, nocheck);
325 //_________________________________________________________________________
326 Int_t AliTRDCalibPadStatus::ProcessEvent(
328 eventHeaderStruct *event,
331 eventHeaderStruct* /*event*/,
338 // process date event
341 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
342 Bool_t result=ProcessEvent(rawReader, nocheck);
346 Fatal("AliTRDCalibPadStatus", "this class was compiled without DATE");
352 //_____________________________________________________________________
353 Bool_t AliTRDCalibPadStatus::TestEvent(Int_t nevent) /*FOLD00*/
357 // fill one oroc and one iroc with random gaus
362 for (Int_t ism=0; ism<18; ism++){
363 for (Int_t ich=0; ich < 5; ich++){
364 for (Int_t ipl=0; ipl < 6; ipl++){
365 for(Int_t irow = 0; irow < fGeo->GetRowMax(ipl,ich,ism); irow++){
366 for(Int_t icol = 0; icol < fGeo->GetColMax(ipl); icol++){
367 for (Int_t iTimeBin=0; iTimeBin<(30*nevent); iTimeBin++){
368 Int_t signal=(Int_t)(ich+8+gRandom->Gaus(0,1.2));
369 if ( signal>0 )Update((ipl+ich*6+ism*6*5),irow,icol,signal,fGeo->GetRowMax(ipl,ich,ism));
379 //_____________________________________________________________________
380 Bool_t AliTRDCalibPadStatus::TestEventHisto(Int_t nevent) /*FOLD00*/
384 // fill one oroc and one iroc with random gaus
389 for (Int_t ism=0; ism<18; ism++){
390 for (Int_t ich=0; ich < 5; ich++){
391 for (Int_t ipl=0; ipl < 6; ipl++){
392 for(Int_t irow = 0; irow < fGeo->GetRowMax(ipl,ich,ism); irow++){
393 for(Int_t icol = 0; icol < fGeo->GetColMax(ipl); icol++){
394 for (Int_t iTimeBin=0; iTimeBin<(30*nevent); iTimeBin++){
395 Int_t signal=(Int_t)(ich+8+gRandom->Gaus(0,1.2));
396 if ( signal>0 )UpdateHisto((ipl+ich*6+ism*6*5),irow,icol,signal,fGeo->GetRowMax(ipl,ich,ism));
406 //_____________________________________________________________________
407 TH2F* AliTRDCalibPadStatus::GetHisto(Int_t det, TObjArray *arr, /*FOLD00*/
408 Int_t nbinsY, Float_t ymin, Float_t ymax,
409 Char_t *type, Bool_t force)
412 // return pointer to histogram
413 // if force is true create a new histogram if it doesn't exist allready
415 if ( !force || arr->UncheckedAt(det) )
416 return (TH2F*)arr->UncheckedAt(det);
418 // if we are forced and histogram doesn't yes exist create it
419 Char_t name[255], title[255];
421 sprintf(name,"hCalib%s%.2d",type,det);
422 sprintf(title,"%s calibration histogram detector %.2d;ADC channel;Channel (pad)",type,det);
425 Int_t nbchannels = fGeo->GetRowMax(GetPlane(det),GetChamber(det),GetSector(det))*fGeo->GetColMax(GetPlane(det));
427 // new histogram with calib information. One value for each pad!
428 TH2F* hist = new TH2F(name,title,
430 nbchannels,0,nbchannels
432 hist->SetDirectory(0);
433 arr->AddAt(hist,det);
437 //_____________________________________________________________________
438 TH2F* AliTRDCalibPadStatus::GetHisto(Int_t det, Bool_t force) /*FOLD00*/
441 // return pointer to histogram
442 // if force is true create a new histogram if it doesn't exist allready
444 TObjArray *arr = &fHistoArray;
445 return GetHisto(det, arr, fAdcMax-fAdcMin, fAdcMin, fAdcMax, "Pedestal", force);
448 //_____________________________________________________________________
449 AliTRDarrayF* AliTRDCalibPadStatus::GetCal(Int_t det, TObjArray* arr, Bool_t force) /*FOLD00*/
452 // return pointer to ROC Calibration
453 // if force is true create a new AliTRDarrayF if it doesn't exist allready
455 if ( !force || arr->UncheckedAt(det) )
456 return (AliTRDarrayF*)arr->UncheckedAt(det);
458 // if we are forced and histogram doesn't yes exist create it
459 AliTRDarrayF *croc = new AliTRDarrayF();
461 Int_t nbpad = fGeo->GetRowMax(GetPlane(det),GetChamber(det),GetSector(det))*fGeo->GetColMax(GetPlane(det));
463 // new AliTRDCalROC. One value for each pad!
465 for(Int_t k = 0; k < nbpad; k++){
468 arr->AddAt(croc,det);
472 //_____________________________________________________________________
473 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRoc(Int_t det, TObjArray* arr, Bool_t force) /*FOLD00*/
476 // return pointer to ROC Calibration
477 // if force is true create a new AliTRDCalROC if it doesn't exist allready
479 if ( !force || arr->UncheckedAt(det) )
480 return (AliTRDCalROC*)arr->UncheckedAt(det);
482 // if we are forced and histogram doesn't yes exist create it
484 // new AliTRDCalROC. One value for each pad!
485 AliTRDCalROC *croc = new AliTRDCalROC(GetPlane(det),GetChamber(det));
486 arr->AddAt(croc,det);
490 //_____________________________________________________________________
491 AliTRDarrayF* AliTRDCalibPadStatus::GetCalEntries(Int_t det, Bool_t force) /*FOLD00*/
494 // return pointer to Carge ROC Calibration
495 // if force is true create a new histogram if it doesn't exist allready
497 TObjArray *arr = &fCalArrayEntries;
498 return GetCal(det, arr, force);
501 //_____________________________________________________________________
502 AliTRDarrayF* AliTRDCalibPadStatus::GetCalMean(Int_t det, Bool_t force) /*FOLD00*/
505 // return pointer to Carge ROC Calibration
506 // if force is true create a new histogram if it doesn't exist allready
508 TObjArray *arr = &fCalArrayMean;
509 return GetCal(det, arr, force);
512 //_____________________________________________________________________
513 AliTRDarrayF* AliTRDCalibPadStatus::GetCalSquares(Int_t det, Bool_t force) /*FOLD00*/
516 // return pointer to Carge ROC Calibration
517 // if force is true create a new histogram if it doesn't exist allready
519 TObjArray *arr = &fCalArraySquares;
520 return GetCal(det, arr, force);
523 //_____________________________________________________________________
524 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocMean(Int_t det, Bool_t force) /*FOLD00*/
527 // return pointer to Carge ROC Calibration
528 // if force is true create a new histogram if it doesn't exist allready
530 TObjArray *arr = &fCalRocArrayMean;
531 return GetCalRoc(det, arr, force);
534 //_____________________________________________________________________
535 AliTRDCalROC* AliTRDCalibPadStatus::GetCalRocRMS(Int_t det, Bool_t force) /*FOLD00*/
538 // return pointer to Carge ROC Calibration
539 // if force is true create a new histogram if it doesn't exist allready
541 TObjArray *arr = &fCalRocArrayRMS;
542 return GetCalRoc(det, arr, force);
545 //_________________________________________________________________________
546 void AliTRDCalibPadStatus::Analyse() /*FOLD00*/
549 // Calcul the rms properly
552 for(Int_t idet = 0; idet < 540; idet++){
555 fCalEntries = ((AliTRDarrayF *)GetCalEntries(idet));
556 fCalMean = ((AliTRDarrayF *)GetCalMean(idet));
557 fCalSquares = ((AliTRDarrayF *)GetCalSquares(idet));
559 if(!fCalEntries) continue;
561 AliTRDCalROC *calRocMean = ((AliTRDCalROC *)GetCalRocMean(idet,kTRUE));
562 AliTRDCalROC *calRocRMS = ((AliTRDCalROC *)GetCalRocRMS(idet,kTRUE));
565 Int_t channels = calRocMean->GetNchannels();
567 for(Int_t ichannels = 0 ; ichannels < channels; ichannels++){
569 Float_t entries = fCalEntries->At(ichannels);
570 Float_t mean = fCalMean->At(ichannels);
571 Float_t squares = fCalSquares->At(ichannels);
575 Double_t rm = TMath::Abs(squares-(mean*mean));
576 rms = TMath::Sqrt(rm);
577 calRocRMS->SetValue(ichannels,rms/10.0);
578 calRocMean->SetValue(ichannels,mean/10.0);
584 //_____________________________________________________________________
585 void AliTRDCalibPadStatus::AnalyseHisto() /*FOLD00*/
588 // Calculate calibration constants
591 Int_t nbinsAdc = fAdcMax-fAdcMin;
599 for (Int_t idet=0; idet<540; idet++){
600 TH2F *hP = GetHisto(idet);
605 AliTRDCalROC *rocMean = GetCalRocMean(idet,kTRUE);
606 AliTRDCalROC *rocRMS = GetCalRocRMS(idet,kTRUE);
608 array_hP = hP->GetArray();
609 Int_t nChannels = rocMean->GetNchannels();
611 for (Int_t iChannel=0; iChannel<nChannels; iChannel++){
612 Int_t offset = (nbinsAdc+2)*(iChannel+1)+1;
613 Double_t ret = AliMathBase::FitGaus(array_hP+offset,nbinsAdc,fAdcMin,fAdcMax,¶m,&dummy);
614 // if the fitting failed set noise and pedestal to 0
615 if ((ret == -4) || (ret == -1)) {
619 if((param[1]/10.0) > 65534.0) param[1] = 0.0;
620 if((param[2]/10.0) > 65534.0) param[2] = 0.0;
621 rocMean->SetValue(iChannel,param[1]/10.0);
622 rocRMS->SetValue(iChannel,param[2]/10.0);
628 //_______________________________________________________________________________________
629 AliTRDCalPadStatus* AliTRDCalibPadStatus::CreateCalPadStatus()
635 AliTRDCalPadStatus* obj = new AliTRDCalPadStatus("padstatus", "padstatus");
637 for (Int_t idet=0; idet<540; ++idet)
639 AliTRDCalSingleChamberStatus *calROC = obj->GetCalROC(idet);
642 fCalEntries = ((AliTRDarrayF *)GetCalEntries(idet));
643 AliTRDCalROC *calRocMean = ((AliTRDCalROC *)GetCalRocMean(idet));
644 AliTRDCalROC *calRocRMS = ((AliTRDCalROC *)GetCalRocRMS(idet));
647 for(Int_t k = 0; k < calROC->GetNchannels(); k++){
648 calROC->SetStatus(k,AliTRDCalPadStatus::kMasked);
650 //printf("no fCalRocMean for %d\n",idet);
655 Int_t channels = calROC->GetNchannels();
659 Float_t meanentries = 0.0;
662 meanentries = GetHisto(idet)->GetEntries()/(channels);
665 else meanentries = TMath::Mean(channels,((TArrayF *)fCalEntries)->GetArray());
666 //Double_t meanmean = calRocMean->GetMean()*10.0;
667 //Double_t meansquares = calRocRMS->GetMean()*10.0;
670 for(Int_t ich = 0; ich < channels; ich++){
672 Float_t entries = 0.0;
675 for(Int_t bin = 0; bin < (fAdcMax-fAdcMin); bin++){
676 entries += GetHisto(idet)->GetArray()[(ich+1)*(fAdcMax-fAdcMin+2)+(bin+1)];
680 else entries = fCalEntries->At(ich);
681 //Float_t mean = calRocMean->GetValue(ich)*10.0;
682 Float_t rms = calRocRMS->GetValue(ich)*10.0;
683 //if(ich > 1720) printf("rms %f\n",rms);
685 //if((entries < 0.3*meanentries) || (TMath::Abs(rms-meansquares) > ((Float_t)(fAdcMax-fAdcMin)/2.0)) || (rms == 0.0)) calROC->SetStatus(ich, AliTRDCalPadStatus::kMasked);
686 if(rms <= 0.0001) calROC->SetStatus(ich, AliTRDCalPadStatus::kMasked);
694 //_____________________________________________________________________
695 void AliTRDCalibPadStatus::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append) /*FOLD00*/
698 // Write class to file
709 TDirectory *backup = gDirectory;
710 TFile f(filename,option.Data());
712 if ( !sDir.IsNull() ){
713 f.mkdir(sDir.Data());
719 if ( backup ) backup->cd();
722 //_____________________________________________________________________________
723 Int_t AliTRDCalibPadStatus::GetPlane(Int_t d) const
726 // Reconstruct the plane number from the detector number
729 return ((Int_t) (d % 6));
733 //_____________________________________________________________________________
734 Int_t AliTRDCalibPadStatus::GetChamber(Int_t d) const
737 // Reconstruct the chamber number from the detector number
740 return ((Int_t) (d % 30) / 6);
744 //_____________________________________________________________________________
745 Int_t AliTRDCalibPadStatus::GetSector(Int_t d) const
748 // Reconstruct the sector number from the detector number
751 return ((Int_t) (d / 30));