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 **************************************************************************/
18 /////////////////////////////////////////////////////////////////////////////////////////
20 // Implementation of the TPC pulser calibration //
22 // Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
24 /////////////////////////////////////////////////////////////////////////////////////////
25 /***************************************************************************
27 ***************************************************************************
29 The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
30 runs performed with the calibration pulser.
32 The information retrieved is
34 - Signal width differences
35 - Amplification variations
37 the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
38 one chip and somewhat large between different chips.
41 For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum]) is created when
42 it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
43 TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
48 Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
49 (see below). These in the end call the Update(...) function.
51 - the Update(...) function:
52 In this function the array fPadSignal is filled with the adc signals between the specified range
53 fFirstTimeBin and fLastTimeBin for the current pad.
54 before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
57 - the ProcessPad() function:
58 Find Pedestal and Noise information
59 - use database information which has to be set by calling
60 SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
61 - if no information from the pedestal data base
62 is available the informaion is calculated on the fly ( see FindPedestal() function )
64 Find the Pulser signal information
65 - calculate mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
66 the Q sum is scaled by pad area
67 (see FindPulserSignal(...) function)
69 Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
70 Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
72 At the end of each event the EndEvent() function is called
74 - the EndEvent() function:
75 calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
76 This is done to overcome syncronisation problems between the trigger and the fec clock.
78 After accumulating the desired statistics the Analyse() function has to be called.
79 - the Analyse() function
80 Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
81 the AliMathBase::GetCOG(...) function, and the calibration
82 storage classes (AliTPCCalROC) are filled for each ROC.
83 The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
88 User interface for filling data:
89 --------------------------------
91 To Fill information one of the following functions can be used:
93 Bool_t ProcessEvent(eventHeaderStruct *event);
95 - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
97 Bool_t ProcessEvent(AliRawReader *rawReader);
98 - process AliRawReader event
99 - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
101 Bool_t ProcessEvent(AliTPCRawStream *rawStream);
102 - process event from AliTPCRawStream
103 - call Update function for signal filling
105 Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
106 iPad, const Int_t iTimeBin, const Float_t signal);
107 - directly fill signal information (sector, row, pad, time bin, pad)
108 to the reference histograms
110 It is also possible to merge two independently taken calibrations using the function
112 void Merge(AliTPCCalibPulser *sig)
113 - copy histograms in 'sig' if the do not exist in this instance
114 - Add histograms in 'sig' to the histograms in this instance if the allready exist
115 - After merging call Analyse again!
119 -- example: filling data using root raw data:
120 void fillSignal(Char_t *filename)
122 rawReader = new AliRawReaderRoot(fileName);
123 if ( !rawReader ) return;
124 AliTPCCalibPulser *calib = new AliTPCCalibPulser;
125 while (rawReader->NextEvent()){
126 calib->ProcessEvent(rawReader);
129 calib->DumpToFile("SignalData.root");
135 What kind of information is stored and how to retrieve them:
136 ------------------------------------------------------------
138 - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
140 TH2F *GetHistoT0(Int_t sector);
141 TH2F *GetHistoRMS(Int_t sector);
142 TH2F *GetHistoQ(Int_t sector);
144 - Accessing the calibration storage objects:
146 AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values
147 AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values
148 AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values
150 example for visualisation:
151 if the file "SignalData.root" was created using the above example one could do the following:
153 TFile fileSignal("SignalData.root")
154 AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
155 sig->GetCalRocT0(0)->Draw("colz");
156 sig->GetCalRocRMS(0)->Draw("colz");
158 or use the AliTPCCalPad functionality:
159 AliTPCCalPad padT0(ped->GetCalPadT0());
160 AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
161 padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
162 padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
167 #include <TObjArray.h>
171 #include <TVectorF.h>
174 #include <TDirectory.h>
180 #include "AliRawReader.h"
181 #include "AliRawReaderRoot.h"
182 #include "AliRawReaderDate.h"
183 #include "AliTPCRawStream.h"
184 #include "AliTPCRawStreamFast.h"
185 #include "AliTPCCalROC.h"
186 #include "AliTPCCalPad.h"
187 #include "AliTPCROC.h"
188 #include "AliTPCParam.h"
189 #include "AliTPCCalibPulser.h"
190 #include "AliTPCcalibDB.h"
191 #include "AliMathBase.h"
193 #include "TTreeStream.h"
201 ClassImp(AliTPCCalibPulser)
203 AliTPCCalibPulser::AliTPCCalibPulser() :
217 fOldRCUformat(kTRUE),
218 fROC(AliTPCROC::Instance()),
220 fParam(new AliTPCParam),
229 fCalRocArrayOutliers(72),
233 fPadTimesArrayEvent(72),
235 fPadRMSArrayEvent(72),
236 fPadPedestalArrayEvent(72),
246 fVTime0OffsetCounter(72),
252 // AliTPCSignal default constructor
256 //_____________________________________________________________________
257 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
259 fFirstTimeBin(sig.fFirstTimeBin),
260 fLastTimeBin(sig.fLastTimeBin),
261 fNbinsT0(sig.fNbinsT0),
262 fXminT0(sig.fXminT0),
263 fXmaxT0(sig.fXmaxT0),
264 fNbinsQ(sig.fNbinsQ),
267 fNbinsRMS(sig.fNbinsRMS),
268 fXminRMS(sig.fXminRMS),
269 fXmaxRMS(sig.fXmaxRMS),
271 fOldRCUformat(kTRUE),
272 fROC(AliTPCROC::Instance()),
274 fParam(new AliTPCParam),
283 fCalRocArrayOutliers(72),
287 fPadTimesArrayEvent(72),
289 fPadRMSArrayEvent(72),
290 fPadPedestalArrayEvent(72),
300 fVTime0OffsetCounter(72),
303 fDebugLevel(sig.fDebugLevel)
306 // AliTPCSignal default constructor
309 for (Int_t iSec = 0; iSec < 72; ++iSec){
310 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
311 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
312 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
313 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
315 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
316 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
317 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
319 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
320 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
321 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
322 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
325 TH2S *hNew = new TH2S(*hQ);
326 hNew->SetDirectory(0);
327 fHistoQArray.AddAt(hNew,iSec);
330 TH2S *hNew = new TH2S(*hT0);
331 hNew->SetDirectory(0);
332 fHistoQArray.AddAt(hNew,iSec);
335 TH2S *hNew = new TH2S(*hRMS);
336 hNew->SetDirectory(0);
337 fHistoQArray.AddAt(hNew,iSec);
342 //_____________________________________________________________________
343 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
346 // assignment operator
348 if (&source == this) return *this;
349 new (this) AliTPCCalibPulser(source);
353 //_____________________________________________________________________
354 AliTPCCalibPulser::~AliTPCCalibPulser()
362 if ( fDebugStreamer) delete fDebugStreamer;
366 void AliTPCCalibPulser::Reset()
369 // Delete all information: Arrays, Histograms, CalRoc objects
371 fCalRocArrayT0.Delete();
372 fCalRocArrayQ.Delete();
373 fCalRocArrayRMS.Delete();
374 fCalRocArrayOutliers.Delete();
376 fHistoQArray.Delete();
377 fHistoT0Array.Delete();
378 fHistoRMSArray.Delete();
380 fPadTimesArrayEvent.Delete();
381 fPadQArrayEvent.Delete();
382 fPadRMSArrayEvent.Delete();
383 fPadPedestalArrayEvent.Delete();
385 //_____________________________________________________________________
386 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
389 const Int_t icTimeBin,
390 const Float_t csignal)
393 // Signal filling methode on the fly pedestal and time offset correction if necessary.
394 // no extra analysis necessary. Assumes knowledge of the signal shape!
395 // assumes that it is looped over consecutive time bins of one pad
398 if (icRow<0) return 0;
399 if (icPad<0) return 0;
400 if (icTimeBin<0) return 0;
401 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
403 if ( icRow<0 || icPad<0 ){
404 AliWarning("Wrong Pad or Row number, skipping!");
408 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
410 //init first pad and sector in this event
411 if ( fCurrentChannel == -1 ) {
412 fCurrentChannel = iChannel;
413 fCurrentSector = icsector;
417 //process last pad if we change to a new one
418 if ( iChannel != fCurrentChannel ){
420 fCurrentChannel = iChannel;
421 fCurrentSector = icsector;
425 //fill signals for current pad
426 fPadSignal[icTimeBin]=csignal;
427 if ( csignal > fMaxPadSignal ){
428 fMaxPadSignal = csignal;
429 fMaxTimeBin = icTimeBin;
433 //_____________________________________________________________________
434 void AliTPCCalibPulser::FindPedestal(Float_t part)
437 // find pedestal and noise for the current pad. Use either database or
438 // truncated mean with part*100%
440 Bool_t noPedestal = kTRUE;;
441 if (fPedestalTPC&&fPadNoiseTPC){
442 //use pedestal database
443 //only load new pedestals if the sector has changed
444 if ( fCurrentSector!=fLastSector ){
445 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
446 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
447 fLastSector=fCurrentSector;
450 if ( fPedestalROC&&fPadNoiseROC ){
451 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
452 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
458 //if we are not running with pedestal database, or for the current sector there is no information
459 //available, calculate the pedestal and noise on the fly
461 const Int_t kPedMax = 100; //maximum pedestal value
470 UShort_t histo[kPedMax];
471 memset(histo,0,kPedMax*sizeof(UShort_t));
473 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
474 padSignal = fPadSignal.GetMatrixArray()[i];
475 if (padSignal<=0) continue;
476 if (padSignal>max && i>10) {
480 if (padSignal>kPedMax-1) continue;
481 histo[Int_t(padSignal+0.5)]++;
485 for (Int_t i=1; i<kPedMax; ++i){
486 if (count1<count0*0.5) median=i;
491 // what if by chance histo[median] == 0 ?!?
492 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
494 for (Int_t idelta=1; idelta<10; ++idelta){
495 if (median-idelta<=0) continue;
496 if (median+idelta>kPedMax) continue;
497 if (count<part*count1){
498 count+=histo[median-idelta];
499 mean +=histo[median-idelta]*(median-idelta);
500 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
501 count+=histo[median+idelta];
502 mean +=histo[median+idelta]*(median+idelta);
503 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
510 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
516 //_____________________________________________________________________
517 void AliTPCCalibPulser::FindPulserSignal(TVectorD ¶m, Float_t &qSum)
520 // Find position, signal width and height of the CE signal (last signal)
521 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
522 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
525 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
526 Int_t cemaxpos = fMaxTimeBin;
527 Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal sum
528 Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal max
529 const Int_t kCemin = 2; // range for the analysis of the ce signal +- channels from the peak
530 const Int_t kCemax = 7;
537 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
538 if ( ceQmax<ceMaxThreshold ) return;
539 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
540 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
541 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
542 ceTime+=signal*(i+0.5);
543 ceRMS +=signal*(i+0.5)*(i+0.5);
548 if (ceQsum>ceSumThreshold) {
550 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
551 //only fill the Time0Offset if pad was not marked as an outlier!
553 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
554 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
556 if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
557 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
558 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
562 //Normalise Q to the pad area
563 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
577 //_____________________________________________________________________
578 void AliTPCCalibPulser::ProcessPad()
581 // Process data of current pad
587 FindPulserSignal(param, qSum);
589 Double_t meanT = param[1];
590 Double_t sigmaT = param[2];
593 //Fill Event T0 counter
594 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
597 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
600 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
603 //Fill debugging info
604 if ( fDebugLevel>0 ){
605 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
606 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
607 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
611 //_____________________________________________________________________
612 void AliTPCCalibPulser::EndEvent()
615 // Process data of current event
618 //check if last pad has allready been processed, if not do so
619 if ( fMaxTimeBin>-1 ) ProcessPad();
621 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
622 for ( Int_t iSec = 0; iSec<72; ++iSec ){
623 TVectorF *vTimes = GetPadTimesEvent(iSec);
624 if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
625 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
626 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
627 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
629 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
630 //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
634 if ( fDebugLevel>0 ){
635 if ( !fDebugStreamer ) {
637 TDirectory *backup = gDirectory;
638 fDebugStreamer = new TTreeSRedirector("deb2.root");
639 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
646 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
647 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
649 UInt_t channel=iChannel;
652 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
653 pad = channel-fROC->GetRowIndexes(sector)[row];
654 padc = pad-(fROC->GetNPads(sector,row)/2);
656 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
657 Form("hSignalD%d.%d.%d",sector,row,pad),
658 fLastTimeBin-fFirstTimeBin,
659 fFirstTimeBin,fLastTimeBin);
662 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
663 h1->Fill(i,fPadSignal(i));
665 (*fDebugStreamer) << "DataPad" <<
666 // "Event=" << fEvent <<
667 "Sector="<< sector <<
671 "PadSec="<< channel <<
685 //_____________________________________________________________________
686 Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
689 // Event Processing loop - AliTPCRawStream
693 Bool_t withInput = kFALSE;
695 while ( rawStreamFast->NextDDL() ){
696 while ( rawStreamFast->NextChannel() ){
697 Int_t isector = rawStreamFast->GetSector(); // current sector
698 Int_t iRow = rawStreamFast->GetRow(); // current row
699 Int_t iPad = rawStreamFast->GetPad(); // current pad
701 while ( rawStreamFast->NextBunch() ){
702 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
703 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
704 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
705 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
706 Update(isector,iRow,iPad,iTimeBin+1,signal);
717 //_____________________________________________________________________
718 Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
721 // Event processing loop - AliRawReader
723 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
724 Bool_t res=ProcessEventFast(rawStreamFast);
725 delete rawStreamFast;
728 //_____________________________________________________________________
729 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
732 // Event Processing loop - AliTPCRawStream
735 rawStream->SetOldRCUFormat(fOldRCUformat);
739 Bool_t withInput = kFALSE;
741 while (rawStream->Next()) {
742 Int_t isector = rawStream->GetSector(); // current sector
743 Int_t iRow = rawStream->GetRow(); // current row
744 Int_t iPad = rawStream->GetPad(); // current pad
745 Int_t iTimeBin = rawStream->GetTime(); // current time bin
746 Float_t signal = rawStream->GetSignal(); // current ADC signal
748 Update(isector,iRow,iPad,iTimeBin,signal);
756 //_____________________________________________________________________
757 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
760 // Event processing loop - AliRawReader
764 AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
766 rawReader->Select("TPC");
768 return ProcessEvent(&rawStream);
770 //_____________________________________________________________________
771 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
774 // Event processing loop - date event
776 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
777 Bool_t result=ProcessEvent(rawReader);
782 //_____________________________________________________________________
783 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
784 Int_t nbinsY, Float_t ymin, Float_t ymax,
785 Char_t *type, Bool_t force)
788 // return pointer to Q histogram
789 // if force is true create a new histogram if it doesn't exist allready
791 if ( !force || arr->UncheckedAt(sector) )
792 return (TH2S*)arr->UncheckedAt(sector);
794 // if we are forced and histogram doesn't yes exist create it
795 Char_t name[255], title[255];
797 sprintf(name,"hCalib%s%.2d",type,sector);
798 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
800 // new histogram with Q calib information. One value for each pad!
801 TH2S* hist = new TH2S(name,title,
803 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
804 hist->SetDirectory(0);
805 arr->AddAt(hist,sector);
808 //_____________________________________________________________________
809 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
812 // return pointer to T0 histogram
813 // if force is true create a new histogram if it doesn't exist allready
815 TObjArray *arr = &fHistoT0Array;
816 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
818 //_____________________________________________________________________
819 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
822 // return pointer to Q histogram
823 // if force is true create a new histogram if it doesn't exist allready
825 TObjArray *arr = &fHistoQArray;
826 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
828 //_____________________________________________________________________
829 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
832 // return pointer to Q histogram
833 // if force is true create a new histogram if it doesn't exist allready
835 TObjArray *arr = &fHistoRMSArray;
836 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
838 //_____________________________________________________________________
839 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
842 // return pointer to Pad Info from 'arr' for the current event and sector
843 // if force is true create it if it doesn't exist allready
845 if ( !force || arr->UncheckedAt(sector) )
846 return (TVectorF*)arr->UncheckedAt(sector);
848 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
849 arr->AddAt(vect,sector);
852 //_____________________________________________________________________
853 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
856 // return pointer to Pad Times Array for the current event and sector
857 // if force is true create it if it doesn't exist allready
859 TObjArray *arr = &fPadTimesArrayEvent;
860 return GetPadInfoEvent(sector,arr,force);
862 //_____________________________________________________________________
863 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
866 // return pointer to Pad Q Array for the current event and sector
867 // if force is true create it if it doesn't exist allready
868 // for debugging purposes only
871 TObjArray *arr = &fPadQArrayEvent;
872 return GetPadInfoEvent(sector,arr,force);
874 //_____________________________________________________________________
875 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
878 // return pointer to Pad RMS Array for the current event and sector
879 // if force is true create it if it doesn't exist allready
880 // for debugging purposes only
882 TObjArray *arr = &fPadRMSArrayEvent;
883 return GetPadInfoEvent(sector,arr,force);
885 //_____________________________________________________________________
886 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
889 // return pointer to Pad RMS Array for the current event and sector
890 // if force is true create it if it doesn't exist allready
891 // for debugging purposes only
893 TObjArray *arr = &fPadPedestalArrayEvent;
894 return GetPadInfoEvent(sector,arr,force);
896 //_____________________________________________________________________
897 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
900 // return pointer to ROC Calibration
901 // if force is true create a new histogram if it doesn't exist allready
903 if ( !force || arr->UncheckedAt(sector) )
904 return (AliTPCCalROC*)arr->UncheckedAt(sector);
906 // if we are forced and histogram doesn't yes exist create it
908 // new AliTPCCalROC for T0 information. One value for each pad!
909 AliTPCCalROC *croc = new AliTPCCalROC(sector);
910 arr->AddAt(croc,sector);
913 //_____________________________________________________________________
914 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
917 // return pointer to Carge ROC Calibration
918 // if force is true create a new histogram if it doesn't exist allready
920 TObjArray *arr = &fCalRocArrayT0;
921 return GetCalRoc(sector, arr, force);
923 //_____________________________________________________________________
924 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
927 // return pointer to T0 ROC Calibration
928 // if force is true create a new histogram if it doesn't exist allready
930 TObjArray *arr = &fCalRocArrayQ;
931 return GetCalRoc(sector, arr, force);
933 //_____________________________________________________________________
934 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
937 // return pointer to signal width ROC Calibration
938 // if force is true create a new histogram if it doesn't exist allready
940 TObjArray *arr = &fCalRocArrayRMS;
941 return GetCalRoc(sector, arr, force);
943 //_____________________________________________________________________
944 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
947 // return pointer to Outliers
948 // if force is true create a new histogram if it doesn't exist allready
950 TObjArray *arr = &fCalRocArrayOutliers;
951 return GetCalRoc(sector, arr, force);
953 //_____________________________________________________________________
954 void AliTPCCalibPulser::ResetEvent()
957 // Reset global counters -- Should be called before each event is processed
966 fPadTimesArrayEvent.Delete();
967 fPadQArrayEvent.Delete();
968 fPadRMSArrayEvent.Delete();
969 fPadPedestalArrayEvent.Delete();
971 for ( Int_t i=0; i<72; ++i ){
973 fVTime0OffsetCounter[i]=0;
976 //_____________________________________________________________________
977 void AliTPCCalibPulser::ResetPad()
980 // Reset pad infos -- Should be called after a pad has been processed
982 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
989 //_____________________________________________________________________
990 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
993 // Merge reference histograms of sig to the current AliTPCCalibPulser
997 for (Int_t iSec=0; iSec<72; ++iSec){
998 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
999 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
1000 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1004 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1005 TH2S *hRefQ = GetHistoQ(iSec);
1006 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1008 TH2S *hist = new TH2S(*hRefQmerge);
1009 hist->SetDirectory(0);
1010 fHistoQArray.AddAt(hist, iSec);
1012 hRefQmerge->SetDirectory(dir);
1015 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1016 TH2S *hRefT0 = GetHistoT0(iSec);
1017 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1019 TH2S *hist = new TH2S(*hRefT0merge);
1020 hist->SetDirectory(0);
1021 fHistoT0Array.AddAt(hist, iSec);
1023 hRefT0merge->SetDirectory(dir);
1025 if ( hRefRMSmerge ){
1026 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1027 TH2S *hRefRMS = GetHistoRMS(iSec);
1028 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1030 TH2S *hist = new TH2S(*hRefRMSmerge);
1031 hist->SetDirectory(0);
1032 fHistoRMSArray.AddAt(hist, iSec);
1034 hRefRMSmerge->SetDirectory(dir);
1039 //_____________________________________________________________________
1040 void AliTPCCalibPulser::Analyse()
1043 // Calculate calibration constants
1047 TVectorD paramT0(3);
1048 TVectorD paramRMS(3);
1049 TMatrixD dummy(3,3);
1051 for (Int_t iSec=0; iSec<72; ++iSec){
1052 TH2S *hT0 = GetHistoT0(iSec);
1053 if (!hT0 ) continue;
1055 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1056 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1057 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1058 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1060 TH2S *hQ = GetHistoQ(iSec);
1061 TH2S *hRMS = GetHistoRMS(iSec);
1063 Short_t *arrayhQ = hQ->GetArray();
1064 Short_t *arrayhT0 = hT0->GetArray();
1065 Short_t *arrayhRMS = hRMS->GetArray();
1067 UInt_t nChannels = fROC->GetNChannels(iSec);
1075 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1077 Float_t cogTime0 = -1000;
1078 Float_t cogQ = -1000;
1079 Float_t cogRMS = -1000;
1082 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1083 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1084 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1086 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,¶mQ,&dummy);
1087 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,¶mT0,&dummy);
1088 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,¶mRMS,&dummy);
1090 cogTime0 = paramT0[1];
1091 cogRMS = paramRMS[1];
1093 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1094 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1095 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1098 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1105 rocQ->SetValue(iChannel, cogQ*cogQ);
1106 rocT0->SetValue(iChannel, cogTime0);
1107 rocRMS->SetValue(iChannel, cogRMS);
1108 rocOut->SetValue(iChannel, cogOut);
1112 if ( fDebugLevel > 0 ){
1113 if ( !fDebugStreamer ) {
1115 TDirectory *backup = gDirectory;
1116 fDebugStreamer = new TTreeSRedirector("deb2.root");
1117 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1120 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1121 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1122 padc = pad-(fROC->GetNPads(iSec,row)/2);
1124 (*fDebugStreamer) << "DataEnd" <<
1125 "Sector=" << iSec <<
1129 "PadSec=" << iChannel <<
1131 "T0=" << cogTime0 <<
1138 delete fDebugStreamer;
1139 fDebugStreamer = 0x0;
1141 //_____________________________________________________________________
1142 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1145 // Write class to file
1154 option = "recreate";
1156 TDirectory *backup = gDirectory;
1157 TFile f(filename,option.Data());
1159 if ( !sDir.IsNull() ){
1160 f.mkdir(sDir.Data());
1166 if ( backup ) backup->cd();
1168 //_____________________________________________________________________
1169 //_________________________ Test Functions ___________________________
1170 //_____________________________________________________________________
1171 TObjArray* AliTPCCalibPulser::TestBinning()
1174 // Function to test the binning of the reference histograms
1175 // type: T0, Q or RMS
1176 // mode: 0 - number of filled bins per channel
1177 // 1 - number of empty bins between filled bins in one ROC
1178 // returns TObjArray with the test histograms type*2+mode:
1179 // position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1182 TObjArray *histArray = new TObjArray(6);
1183 const Char_t *type[] = {"T0","Q","RMS"};
1184 Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1186 for (Int_t itype = 0; itype<3; ++itype){
1187 for (Int_t imode=0; imode<2; ++imode){
1188 Int_t icount = itype*2+imode;
1189 histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1190 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1199 for (Int_t itype = 0; itype<3; ++itype){
1200 for (Int_t iSec=0; iSec<72; ++iSec){
1201 if ( itype == 0 ) hRef = GetHistoT0(iSec);
1202 if ( itype == 1 ) hRef = GetHistoQ(iSec);
1203 if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1204 if ( hRef == 0x0 ) continue;
1205 array = (hRef->GetArray());
1206 UInt_t nChannels = fROC->GetNChannels(iSec);
1209 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1211 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1214 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1215 if ( array[offset+iBin]>0 ) {
1217 if ( c1 && c2 ) nempty++;
1220 else if ( c1 ) c2 = 1;
1223 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1225 ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);