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 **************************************************************************/
20 //-------------------------------------------------------
21 // Implementation of the TPC pulser calibration
23 // Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
26 //-------------------------------------------------------
34 #include <TObjArray.h>
46 #include <TStopwatch.h>
49 #include <TDirectory.h>
54 #include "AliRawReader.h"
55 #include "AliRawReaderRoot.h"
56 #include "AliTPCRawStream.h"
57 #include "AliTPCCalROC.h"
58 #include "AliTPCCalPad.h"
59 #include "AliTPCROC.h"
60 #include "AliTPCParam.h"
61 #include "AliTPCCalibSignal.h"
62 #include "AliTPCcalibDB.h"
64 #include "TTreeStream.h"
65 ClassImp(AliTPCCalibSignal) /*FOLD00*/
67 AliTPCCalibSignal::AliTPCCalibSignal() : /*FOLD00*/
74 fROC(AliTPCROC::Instance()),
75 fParam(new AliTPCParam),
81 fCalRocArrayOutliers(72),
85 fPadTimesArrayEvent(72),
87 fPadRMSArrayEvent(72),
88 fPadPedestalArrayEvent(72),
96 fVTime0Offset1Counter(72),
102 // AliTPCSignal default constructor
106 TDirectory *backup = gDirectory;
107 fDebugStreamer = new TTreeSRedirector("deb2.root");
108 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
111 //_____________________________________________________________________
112 AliTPCCalibSignal::AliTPCCalibSignal(const AliTPCCalibSignal &sig) :
114 fFirstTimeBin(sig.fFirstTimeBin),
115 fLastTimeBin(sig.fLastTimeBin),
116 fFirstTimeBinT0(sig.fFirstTimeBinT0),
117 fLastTimeBinT0(sig.fLastTimeBinT0),
119 fROC(AliTPCROC::Instance()),
120 fParam(new AliTPCParam),
121 fPedestalTPC(sig.fPedestalTPC),
122 fBpedestal(sig.fBpedestal),
126 fCalRocArrayOutliers(72),
130 fPadTimesArrayEvent(72),
132 fPadRMSArrayEvent(72),
133 fPadPedestalArrayEvent(72),
141 fVTime0Offset1Counter(72),
144 fDebugLevel(sig.fDebugLevel)
147 // AliTPCSignal default constructor
150 for (Int_t iSec = 0; iSec < 72; iSec++){
151 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
152 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
153 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
154 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
156 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
157 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
158 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
160 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
161 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
162 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
163 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
166 TH2S *hNew = new TH2S(*hQ);
167 hNew->SetDirectory(0);
168 fHistoQArray.AddAt(hNew,iSec);
171 TH2S *hNew = new TH2S(*hT0);
172 hNew->SetDirectory(0);
173 fHistoQArray.AddAt(hNew,iSec);
176 TH2S *hNew = new TH2S(*hRMS);
177 hNew->SetDirectory(0);
178 fHistoQArray.AddAt(hNew,iSec);
184 TDirectory *backup = gDirectory;
185 fDebugStreamer = new TTreeSRedirector("deb2.root");
186 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
189 //_____________________________________________________________________
190 AliTPCCalibSignal& AliTPCCalibSignal::operator = (const AliTPCCalibSignal &source)
193 // assignment operator
195 if (&source == this) return *this;
196 new (this) AliTPCCalibSignal(source);
200 //_____________________________________________________________________
201 AliTPCCalibSignal::~AliTPCCalibSignal()
207 fCalRocArrayT0.Delete();
208 fCalRocArrayQ.Delete();
209 fCalRocArrayRMS.Delete();
211 fHistoQArray.Delete();
212 fHistoT0Array.Delete();
213 fHistoRMSArray.Delete();
215 fPadTimesArrayEvent.Delete();
216 fPadQArrayEvent.Delete();
217 fPadRMSArrayEvent.Delete();
218 fPadPedestalArrayEvent.Delete();
220 if ( fDebugStreamer) delete fDebugStreamer;
224 //_____________________________________________________________________
225 Int_t AliTPCCalibSignal::Update(const Int_t icsector, /*FOLD00*/
228 const Int_t icTimeBin,
229 const Float_t csignal)
232 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
233 // no extra analysis necessary. Assumes knowledge of the signal shape!
234 // assumes that it is looped over consecutive time bins of one pad
236 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
238 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
240 //init first pad and sector in this event
241 if ( fCurrentChannel == -1 ) {
242 fCurrentChannel = iChannel;
243 fCurrentSector = icsector;
247 //process last pad if we change to a new one
248 if ( iChannel != fCurrentChannel ){
250 fCurrentChannel = iChannel;
251 fCurrentSector = icsector;
255 if ( fDebugLevel > 6 )
256 if ( icTimeBin == fLastTimeBin )
257 printf("Filling sector: %d, channel: %d\n",icsector,iChannel);
259 //fill signals for current pad
260 fPadSignal[icTimeBin]=csignal;
261 if ( csignal > fMaxPadSignal ){
262 fMaxPadSignal = csignal;
263 fMaxTimeBin = icTimeBin;
267 //_____________________________________________________________________
268 void AliTPCCalibSignal::ProcessPad() /*FOLD00*/
271 // Process data of current pad
273 if ( fDebugLevel > 4 )
274 printf("process: sector-pad: %.2d-%.4d ... ", fCurrentSector, fCurrentChannel);
276 Float_t pedestal = 0;
279 //!!!!!!! does not work like this
280 //use pedestal database
281 AliTPCCalROC *pedestalROC = 0x0;
283 //only load new pedestals if the sector has changed
284 if ( fCurrentSector!=fLastSector ){
285 pedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
286 fLastSector=fCurrentSector;
289 pedestal = pedestalROC->GetValue(fCurrentChannel);
292 if ( fDebugLevel > 4 )
293 printf("pedestals ... ");
295 //find pedestal for pad on the fly
296 //using a few timebins before the signal
297 Int_t pminus1 = 10, pminus2=5;
300 for (Int_t i=fMaxTimeBin-pminus1; i<fMaxTimeBin-pminus2+1; i++){
301 if ( i>fFirstTimeBin && i<fLastTimeBin ){
302 pedestal+=fPadSignal[i];
307 if ( sumN>0 ) pedestal/=sumN;
308 if ( fDebugLevel > 4 )
309 printf("%.2f ... ",pedestal);
315 //find signal mean and sigma
316 Int_t tminus = 2, tplus=7;
317 Double_t meanT=0, sigmaT=0, Qsum=0;
319 if ( fDebugLevel > 4 )
320 printf(" mean +- sigma (sum) ... ");
322 for (Int_t i=fMaxTimeBin-tminus; i<fMaxTimeBin+tplus; i++){
323 if ( i>fFirstTimeBin && i<fLastTimeBin ){
324 Double_t val=fPadSignal[i]-pedestal;
325 meanT+=val*(i+.5); //+.5: center of the timebin
326 sigmaT+=val*(i+.5)*(i+.5);
333 //!!!! What to do if Qsum == 0???
334 //!!!! Should there be some threshold for max - pedestal and/or Qsum???
335 //!!!! What if Qsum < 0
336 //!!!! only fill time0 offset if Qsum > 0???
340 sigmaT = TMath::Sqrt(TMath::Abs(meanT*meanT - sigmaT));
342 //fill Time0 offset data for this event
343 fVTime0Offset1[fCurrentSector]+=meanT;
344 fVTime0Offset1Counter[fCurrentSector]++;
347 meanT = fLastTimeBinT0+1; //put to overflow bin
348 sigmaT = fLastTimeBinT0-fFirstTimeBinT0; //put to overflow bin
351 if ( fDebugLevel > 4 )
352 printf("%.3f +- %.3f (%.3f) ... ",meanT,sigmaT,Qsum);
355 if ( fDebugLevel > 4 )
356 printf("filling ... ");
359 //Fill Event T0 counter
360 (*GetPadTimesEvent(fCurrentSector,kTRUE))[fCurrentChannel] = meanT;
363 //Normalise Q to pad area of irocs
364 Float_t norm = fParam->GetPadPitchWidth(0)*fParam->GetPadPitchLength(0,0)/(
365 fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow));
367 // if (fCurrentChannel == 0) printf("sec, norm: %d, %f\n", fCurrentSector, norm);
370 GetHistoQ(fCurrentSector,kTRUE)->Fill(fCurrentChannel, TMath::Sqrt(Qsum*norm));
373 GetHistoRMS(fCurrentSector,kTRUE)->Fill(fCurrentChannel, sigmaT);
376 //Fill debugging info
377 if ( fDebugLevel>0 ){
378 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=pedestal;
379 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
380 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=Qsum;
383 if ( fDebugLevel > 4 )
384 printf("reset pad ...");
388 if ( fDebugLevel > 4 )
391 //_____________________________________________________________________
392 void AliTPCCalibSignal::EndEvent() /*FOLD00*/
395 // Process data of current pad
397 printf("end event\n");
398 for ( Int_t iSec = 0; iSec<72; iSec++ ){
399 TVectorF *vTimes = GetPadTimesEvent(iSec);
400 if ( !vTimes ) continue;
402 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++ ){
403 Float_t Time0 = fVTime0Offset1[iSec]/fVTime0Offset1Counter[iSec];
404 Float_t Time = (*vTimes)[iChannel];
406 GetHistoT0(iSec,kTRUE)->Fill(iChannel,Time-Time0);
410 if ( fDebugLevel>0 ){
415 Float_t Q = (*GetPadQEvent(iSec))[iChannel];
416 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
418 UInt_t channel=iChannel;
421 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
422 pad = channel-fROC->GetRowIndexes(sector)[row];
423 padc = pad-(fROC->GetNPads(sector,row)/2);
425 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
426 Form("hSignalD%d.%d.%d",sector,row,pad),
427 fLastTimeBin-fFirstTimeBin,
428 fFirstTimeBin,fLastTimeBin);
431 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
432 h1->Fill(i,fPadSignal(i));
434 (*fDebugStreamer) << "DataPad" <<
435 "Event=" << fEvent <<
436 "Sector="<< sector <<
440 "PadSec="<< channel <<
456 //_____________________________________________________________________
457 Bool_t AliTPCCalibSignal::ProcessEvent(AliRawReader *rawReader) /*FOLD00*/
464 AliTPCRawStream input(rawReader);
466 rawReader->Select("TPC");
468 input.SetOldRCUFormat(1);
469 printf("start event processing\n");
473 Bool_t withInput = kFALSE;
475 while (input.Next()) {
477 Int_t isector = input.GetSector(); // current sector
478 Int_t iRow = input.GetRow(); // current row
479 Int_t iPad = input.GetPad(); // current pad
480 Int_t iTimeBin = input.GetTime(); // current time bin
481 Float_t signal = input.GetSignal(); // current ADC signal
483 Update(isector,iRow,iPad,iTimeBin,signal);
492 printf("end event processing\n");
494 fDebugStreamer->GetFile()->Write();
497 //_____________________________________________________________________
498 TH2S* AliTPCCalibSignal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
499 Int_t nbinsY, Float_t ymin, Float_t ymax,
500 Char_t *type, Bool_t force)
503 // return pointer to Q histogram
504 // if force is true create a new histogram if it doesn't exist allready
506 if ( !force || arr->UncheckedAt(sector) )
507 return (TH2S*)arr->UncheckedAt(sector);
509 // if we are forced and histogram doesn't yes exist create it
510 Char_t name[255], title[255];
512 sprintf(name,"hCalib%s%.2d",type,sector);
513 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
515 // new histogram with Q calib information. One value for each pad!
516 TH2S* hist = new TH2S(name,title,
517 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector),
519 hist->SetDirectory(0);
520 arr->AddAt(hist,sector);
523 //_____________________________________________________________________
524 TH2S* AliTPCCalibSignal::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
527 // return pointer to T0 histogram
528 // if force is true create a new histogram if it doesn't exist allready
530 TObjArray *arr = &fHistoT0Array;
531 return GetHisto(sector, arr, 200, -2, 2, "T0", force);
533 //_____________________________________________________________________
534 TH2S* AliTPCCalibSignal::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
537 // return pointer to Q histogram
538 // if force is true create a new histogram if it doesn't exist allready
540 TObjArray *arr = &fHistoQArray;
541 return GetHisto(sector, arr, 150, 24, 55, "Q", force);
543 //_____________________________________________________________________
544 TH2S* AliTPCCalibSignal::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
547 // return pointer to Q histogram
548 // if force is true create a new histogram if it doesn't exist allready
550 TObjArray *arr = &fHistoRMSArray;
551 return GetHisto(sector, arr, 100, 0, 5, "RMS", force);
553 //_____________________________________________________________________
554 TVectorF* AliTPCCalibSignal::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
557 // return pointer to Pad Info from 'arr' for the current event and sector
558 // if force is true create it if it doesn't exist allready
560 if ( !force || arr->UncheckedAt(sector) )
561 return (TVectorF*)arr->UncheckedAt(sector);
563 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
564 arr->AddAt(vect,sector);
567 //_____________________________________________________________________
568 TVectorF* AliTPCCalibSignal::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
571 // return pointer to Pad Times Array for the current event and sector
572 // if force is true create it if it doesn't exist allready
574 TObjArray *arr = &fPadTimesArrayEvent;
575 return GetPadInfoEvent(sector,arr,force);
577 //_____________________________________________________________________
578 TVectorF* AliTPCCalibSignal::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
581 // return pointer to Pad Q Array for the current event and sector
582 // if force is true create it if it doesn't exist allready
583 // for debugging purposes only
586 TObjArray *arr = &fPadQArrayEvent;
587 return GetPadInfoEvent(sector,arr,force);
589 //_____________________________________________________________________
590 TVectorF* AliTPCCalibSignal::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
593 // return pointer to Pad RMS Array for the current event and sector
594 // if force is true create it if it doesn't exist allready
595 // for debugging purposes only
597 TObjArray *arr = &fPadRMSArrayEvent;
598 return GetPadInfoEvent(sector,arr,force);
600 //_____________________________________________________________________
601 TVectorF* AliTPCCalibSignal::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
604 // return pointer to Pad RMS Array for the current event and sector
605 // if force is true create it if it doesn't exist allready
606 // for debugging purposes only
608 TObjArray *arr = &fPadPedestalArrayEvent;
609 return GetPadInfoEvent(sector,arr,force);
611 //_____________________________________________________________________
612 AliTPCCalROC* AliTPCCalibSignal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
615 // return pointer to ROC Calibration
616 // if force is true create a new histogram if it doesn't exist allready
618 if ( !force || arr->UncheckedAt(sector) )
619 return (AliTPCCalROC*)arr->UncheckedAt(sector);
621 // if we are forced and histogram doesn't yes exist create it
623 // new AliTPCCalROC for T0 information. One value for each pad!
624 AliTPCCalROC *croc = new AliTPCCalROC(sector);
626 for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
627 croc->SetValue(iChannel, 0);
629 arr->AddAt(croc,sector);
632 //_____________________________________________________________________
633 AliTPCCalROC* AliTPCCalibSignal::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
636 // return pointer to Carge ROC Calibration
637 // if force is true create a new histogram if it doesn't exist allready
639 TObjArray *arr = &fCalRocArrayT0;
640 return GetCalRoc(sector, arr, force);
642 //_____________________________________________________________________
643 AliTPCCalROC* AliTPCCalibSignal::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
646 // return pointer to T0 ROC Calibration
647 // if force is true create a new histogram if it doesn't exist allready
649 TObjArray *arr = &fCalRocArrayQ;
650 return GetCalRoc(sector, arr, force);
652 //_____________________________________________________________________
653 AliTPCCalROC* AliTPCCalibSignal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
656 // return pointer to signal width ROC Calibration
657 // if force is true create a new histogram if it doesn't exist allready
659 TObjArray *arr = &fCalRocArrayRMS;
660 return GetCalRoc(sector, arr, force);
662 //_____________________________________________________________________
663 AliTPCCalROC* AliTPCCalibSignal::GetCalRocOutliers(Int_t sector, Bool_t force)
666 // return pointer to Outliers
667 // if force is true create a new histogram if it doesn't exist allready
669 TObjArray *arr = &fCalRocArrayOutliers;
670 return GetCalRoc(sector, arr, force);
672 //_____________________________________________________________________
673 void AliTPCCalibSignal::ResetEvent() /*FOLD00*/
676 // Reset global counters -- Should be called before each event is processed
685 fPadTimesArrayEvent.Delete();
686 fPadQArrayEvent.Delete();
687 fPadRMSArrayEvent.Delete();
688 fPadPedestalArrayEvent.Delete();
690 for ( Int_t i=0; i<72; i++ ){
692 fVTime0Offset1Counter[i]=0;
695 //_____________________________________________________________________
696 void AliTPCCalibSignal::ResetPad() /*FOLD00*/
699 // Reset pad infos -- Should be called after a pad has been processed
701 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
706 //_____________________________________________________________________
707 void AliTPCCalibSignal::Analyse() /*FOLD00*/
710 // Calculate calibration constants
714 for (Int_t iSec=0; iSec<72; iSec++){
715 TH2S *hT0 = GetHistoT0(iSec);
718 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
719 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
720 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
721 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
723 TH2S *hQ = GetHistoQ(iSec);
724 TH2S *hRMS = GetHistoRMS(iSec);
732 for (UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++){
735 Float_t cogTime0 = -1000;
736 Float_t cogQ = -1000;
737 Float_t cogRMS = -1000;
740 hQ->SetAxisRange(iChannel,iChannel);
741 hT0->SetAxisRange(iChannel,iChannel);
742 hRMS->SetAxisRange(iChannel,iChannel);
744 cogTime0 = hT0->GetMean(2);
745 cogQ = hQ->GetMean(2);
746 cogRMS = hRMS->GetMean(2);
748 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
755 rocQ->SetValue(iChannel, cogQ*cogQ);
756 rocT0->SetValue(iChannel, cogTime0);
757 rocRMS->SetValue(iChannel, cogRMS);
758 rocOut->SetValue(iChannel, cogOut);
762 if ( fDebugLevel > 0 ){
763 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
764 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
765 padc = pad-(fROC->GetNPads(iSec,row)/2);
767 (*fDebugStreamer) << "DataEnd" <<
772 "PadSec=" << iChannel <<
783 delete fDebugStreamer;
784 fDebugStreamer = 0x0;
786 //_____________________________________________________________________
787 void AliTPCCalibSignal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
790 // Write class to file
793 TDirectory *backup = gDirectory;
802 TFile f(filename,option.Data());
803 if ( !sDir.IsNull() ){
804 f.mkdir(sDir.Data());
807 gDirectory->WriteTObject(this);
810 if ( backup ) backup->cd();