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() :
216 fIsZeroSuppressed(kFALSE),
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
257 //_____________________________________________________________________
258 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
260 fFirstTimeBin(sig.fFirstTimeBin),
261 fLastTimeBin(sig.fLastTimeBin),
262 fNbinsT0(sig.fNbinsT0),
263 fXminT0(sig.fXminT0),
264 fXmaxT0(sig.fXmaxT0),
265 fNbinsQ(sig.fNbinsQ),
268 fNbinsRMS(sig.fNbinsRMS),
269 fXminRMS(sig.fXminRMS),
270 fXmaxRMS(sig.fXmaxRMS),
271 fIsZeroSuppressed(sig.fIsZeroSuppressed),
273 fROC(AliTPCROC::Instance()),
275 fParam(new AliTPCParam),
284 fCalRocArrayOutliers(72),
288 fPadTimesArrayEvent(72),
290 fPadRMSArrayEvent(72),
291 fPadPedestalArrayEvent(72),
301 fVTime0OffsetCounter(72),
304 fDebugLevel(sig.fDebugLevel)
307 // AliTPCSignal default constructor
310 for (Int_t iSec = 0; iSec < 72; ++iSec){
311 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
312 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
313 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
314 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
316 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
317 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
318 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
320 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
321 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
322 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
323 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
326 TH2S *hNew = new TH2S(*hQ);
327 hNew->SetDirectory(0);
328 fHistoQArray.AddAt(hNew,iSec);
331 TH2S *hNew = new TH2S(*hT0);
332 hNew->SetDirectory(0);
333 fHistoQArray.AddAt(hNew,iSec);
336 TH2S *hNew = new TH2S(*hRMS);
337 hNew->SetDirectory(0);
338 fHistoQArray.AddAt(hNew,iSec);
344 //_____________________________________________________________________
345 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
348 // assignment operator
350 if (&source == this) return *this;
351 new (this) AliTPCCalibPulser(source);
355 //_____________________________________________________________________
356 AliTPCCalibPulser::~AliTPCCalibPulser()
364 if ( fDebugStreamer) delete fDebugStreamer;
368 void AliTPCCalibPulser::Reset()
371 // Delete all information: Arrays, Histograms, CalRoc objects
373 fCalRocArrayT0.Delete();
374 fCalRocArrayQ.Delete();
375 fCalRocArrayRMS.Delete();
376 fCalRocArrayOutliers.Delete();
378 fHistoQArray.Delete();
379 fHistoT0Array.Delete();
380 fHistoRMSArray.Delete();
382 fPadTimesArrayEvent.Delete();
383 fPadQArrayEvent.Delete();
384 fPadRMSArrayEvent.Delete();
385 fPadPedestalArrayEvent.Delete();
387 //_____________________________________________________________________
388 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
391 const Int_t icTimeBin,
392 const Float_t csignal)
395 // Signal filling methode on the fly pedestal and time offset correction if necessary.
396 // no extra analysis necessary. Assumes knowledge of the signal shape!
397 // assumes that it is looped over consecutive time bins of one pad
400 if (icRow<0) return 0;
401 if (icPad<0) return 0;
402 if (icTimeBin<0) return 0;
403 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
405 if ( icRow<0 || icPad<0 ){
406 AliWarning("Wrong Pad or Row number, skipping!");
410 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
412 //init first pad and sector in this event
413 if ( fCurrentChannel == -1 ) {
414 fCurrentChannel = iChannel;
415 fCurrentSector = icsector;
419 //process last pad if we change to a new one
420 if ( iChannel != fCurrentChannel ){
422 fCurrentChannel = iChannel;
423 fCurrentSector = icsector;
427 //fill signals for current pad
428 fPadSignal[icTimeBin]=csignal;
429 if ( csignal > fMaxPadSignal ){
430 fMaxPadSignal = csignal;
431 fMaxTimeBin = icTimeBin;
435 //_____________________________________________________________________
436 void AliTPCCalibPulser::FindPedestal(Float_t part)
439 // find pedestal and noise for the current pad. Use either database or
440 // truncated mean with part*100%
442 Bool_t noPedestal = kTRUE;;
443 if (fPedestalTPC&&fPadNoiseTPC){
444 //use pedestal database
445 //only load new pedestals if the sector has changed
446 if ( fCurrentSector!=fLastSector ){
447 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
448 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
449 fLastSector=fCurrentSector;
452 if ( fPedestalROC&&fPadNoiseROC ){
453 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
454 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
460 //if we are not running with pedestal database, or for the current sector there is no information
461 //available, calculate the pedestal and noise on the fly
463 const Int_t kPedMax = 100; //maximum pedestal value
472 UShort_t histo[kPedMax];
473 memset(histo,0,kPedMax*sizeof(UShort_t));
475 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
476 padSignal = fPadSignal.GetMatrixArray()[i];
477 if (padSignal<=0) continue;
478 if (padSignal>max && i>10) {
482 if (padSignal>kPedMax-1) continue;
483 histo[Int_t(padSignal+0.5)]++;
487 for (Int_t i=1; i<kPedMax; ++i){
488 if (count1<count0*0.5) median=i;
493 // what if by chance histo[median] == 0 ?!?
494 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
496 for (Int_t idelta=1; idelta<10; ++idelta){
497 if (median-idelta<=0) continue;
498 if (median+idelta>kPedMax) continue;
499 if (count<part*count1){
500 count+=histo[median-idelta];
501 mean +=histo[median-idelta]*(median-idelta);
502 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
503 count+=histo[median+idelta];
504 mean +=histo[median+idelta]*(median+idelta);
505 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
512 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
517 fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
519 //_____________________________________________________________________
520 void AliTPCCalibPulser::FindPulserSignal(TVectorD ¶m, Float_t &qSum)
523 // Find position, signal width and height of the CE signal (last signal)
524 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
525 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
528 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
529 Int_t cemaxpos = fMaxTimeBin;
530 Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal sum
531 Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal max
532 const Int_t kCemin = 2; // range for the analysis of the ce signal +- channels from the peak
533 const Int_t kCemax = 7;
540 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
541 if ( ceQmax<ceMaxThreshold ) return;
542 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
543 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
544 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
545 ceTime+=signal*(i+0.5);
546 ceRMS +=signal*(i+0.5)*(i+0.5);
551 if (ceQsum>ceSumThreshold) {
553 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
554 //only fill the Time0Offset if pad was not marked as an outlier!
556 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
557 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
559 if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
560 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
561 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
565 //Normalise Q to the pad area
566 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
580 //_____________________________________________________________________
581 void AliTPCCalibPulser::ProcessPad()
584 // Process data of current pad
590 FindPulserSignal(param, qSum);
592 Double_t meanT = param[1];
593 Double_t sigmaT = param[2];
596 //Fill Event T0 counter
597 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
600 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
603 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
606 //Fill debugging info
607 if ( fDebugLevel>0 ){
608 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
609 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
610 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
614 //_____________________________________________________________________
615 void AliTPCCalibPulser::EndEvent()
618 // Process data of current event
621 //check if last pad has allready been processed, if not do so
622 if ( fMaxTimeBin>-1 ) ProcessPad();
624 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
625 for ( Int_t iSec = 0; iSec<72; ++iSec ){
626 TVectorF *vTimes = GetPadTimesEvent(iSec);
627 if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
628 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
629 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
630 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
632 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
633 //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
637 if ( fDebugLevel>0 ){
638 if ( !fDebugStreamer ) {
640 TDirectory *backup = gDirectory;
641 fDebugStreamer = new TTreeSRedirector("deb2.root");
642 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
649 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
650 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
652 UInt_t channel=iChannel;
655 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
656 pad = channel-fROC->GetRowIndexes(sector)[row];
657 padc = pad-(fROC->GetNPads(sector,row)/2);
659 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
660 Form("hSignalD%d.%d.%d",sector,row,pad),
661 fLastTimeBin-fFirstTimeBin,
662 fFirstTimeBin,fLastTimeBin);
665 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
666 h1->Fill(i,fPadSignal(i));
668 (*fDebugStreamer) << "DataPad" <<
669 // "Event=" << fEvent <<
670 "Sector="<< sector <<
674 "PadSec="<< channel <<
688 //_____________________________________________________________________
689 Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
692 // Event Processing loop - AliTPCRawStream
695 Bool_t withInput = kFALSE;
696 while ( rawStreamFast->NextDDL() ){
697 while ( rawStreamFast->NextChannel() ){
698 Int_t isector = rawStreamFast->GetSector(); // current sector
699 Int_t iRow = rawStreamFast->GetRow(); // current row
700 Int_t iPad = rawStreamFast->GetPad(); // current pad
702 while ( rawStreamFast->NextBunch() ){
703 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
704 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
705 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
706 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
707 Update(isector,iRow,iPad,iTimeBin+1,signal);
718 //_____________________________________________________________________
719 Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
722 // Event processing loop - AliRawReader
724 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
725 Bool_t res=ProcessEventFast(rawStreamFast);
726 delete rawStreamFast;
729 //_____________________________________________________________________
730 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
733 // Event Processing loop - AliTPCRawStream
738 Bool_t withInput = kFALSE;
740 while (rawStream->Next()) {
741 Int_t isector = rawStream->GetSector(); // current sector
742 Int_t iRow = rawStream->GetRow(); // current row
743 Int_t iPad = rawStream->GetPad(); // current pad
744 Int_t iTimeBin = rawStream->GetTime(); // current time bin
745 Float_t signal = rawStream->GetSignal(); // current ADC signal
747 Update(isector,iRow,iPad,iTimeBin,signal);
755 //_____________________________________________________________________
756 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
759 // Event processing loop - AliRawReader
763 AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
765 rawReader->Select("TPC");
767 return ProcessEvent(&rawStream);
769 //_____________________________________________________________________
770 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
773 // Event processing loop - date event
775 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
776 Bool_t result=ProcessEvent(rawReader);
781 //_____________________________________________________________________
782 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
783 Int_t nbinsY, Float_t ymin, Float_t ymax,
784 Char_t *type, Bool_t force)
787 // return pointer to Q histogram
788 // if force is true create a new histogram if it doesn't exist allready
790 if ( !force || arr->UncheckedAt(sector) )
791 return (TH2S*)arr->UncheckedAt(sector);
793 // if we are forced and histogram doesn't yes exist create it
794 Char_t name[255], title[255];
796 sprintf(name,"hCalib%s%.2d",type,sector);
797 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
799 // new histogram with Q calib information. One value for each pad!
800 TH2S* hist = new TH2S(name,title,
802 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
803 hist->SetDirectory(0);
804 arr->AddAt(hist,sector);
807 //_____________________________________________________________________
808 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
811 // return pointer to T0 histogram
812 // if force is true create a new histogram if it doesn't exist allready
814 TObjArray *arr = &fHistoT0Array;
815 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
817 //_____________________________________________________________________
818 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
821 // return pointer to Q histogram
822 // if force is true create a new histogram if it doesn't exist allready
824 TObjArray *arr = &fHistoQArray;
825 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
827 //_____________________________________________________________________
828 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
831 // return pointer to Q histogram
832 // if force is true create a new histogram if it doesn't exist allready
834 TObjArray *arr = &fHistoRMSArray;
835 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
837 //_____________________________________________________________________
838 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
841 // return pointer to Pad Info from 'arr' for the current event and sector
842 // if force is true create it if it doesn't exist allready
844 if ( !force || arr->UncheckedAt(sector) )
845 return (TVectorF*)arr->UncheckedAt(sector);
847 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
848 arr->AddAt(vect,sector);
851 //_____________________________________________________________________
852 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
855 // return pointer to Pad Times Array for the current event and sector
856 // if force is true create it if it doesn't exist allready
858 TObjArray *arr = &fPadTimesArrayEvent;
859 return GetPadInfoEvent(sector,arr,force);
861 //_____________________________________________________________________
862 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
865 // return pointer to Pad Q Array for the current event and sector
866 // if force is true create it if it doesn't exist allready
867 // for debugging purposes only
870 TObjArray *arr = &fPadQArrayEvent;
871 return GetPadInfoEvent(sector,arr,force);
873 //_____________________________________________________________________
874 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
877 // return pointer to Pad RMS Array for the current event and sector
878 // if force is true create it if it doesn't exist allready
879 // for debugging purposes only
881 TObjArray *arr = &fPadRMSArrayEvent;
882 return GetPadInfoEvent(sector,arr,force);
884 //_____________________________________________________________________
885 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
888 // return pointer to Pad RMS Array for the current event and sector
889 // if force is true create it if it doesn't exist allready
890 // for debugging purposes only
892 TObjArray *arr = &fPadPedestalArrayEvent;
893 return GetPadInfoEvent(sector,arr,force);
895 //_____________________________________________________________________
896 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
899 // return pointer to ROC Calibration
900 // if force is true create a new histogram if it doesn't exist allready
902 if ( !force || arr->UncheckedAt(sector) )
903 return (AliTPCCalROC*)arr->UncheckedAt(sector);
905 // if we are forced and histogram doesn't yes exist create it
907 // new AliTPCCalROC for T0 information. One value for each pad!
908 AliTPCCalROC *croc = new AliTPCCalROC(sector);
909 arr->AddAt(croc,sector);
912 //_____________________________________________________________________
913 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
916 // return pointer to Carge ROC Calibration
917 // if force is true create a new histogram if it doesn't exist allready
919 TObjArray *arr = &fCalRocArrayT0;
920 return GetCalRoc(sector, arr, force);
922 //_____________________________________________________________________
923 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
926 // return pointer to T0 ROC Calibration
927 // if force is true create a new histogram if it doesn't exist allready
929 TObjArray *arr = &fCalRocArrayQ;
930 return GetCalRoc(sector, arr, force);
932 //_____________________________________________________________________
933 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
936 // return pointer to signal width ROC Calibration
937 // if force is true create a new histogram if it doesn't exist allready
939 TObjArray *arr = &fCalRocArrayRMS;
940 return GetCalRoc(sector, arr, force);
942 //_____________________________________________________________________
943 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
946 // return pointer to Outliers
947 // if force is true create a new histogram if it doesn't exist allready
949 TObjArray *arr = &fCalRocArrayOutliers;
950 return GetCalRoc(sector, arr, force);
952 //_____________________________________________________________________
953 void AliTPCCalibPulser::ResetEvent()
956 // Reset global counters -- Should be called before each event is processed
965 fPadTimesArrayEvent.Delete();
966 fPadQArrayEvent.Delete();
967 fPadRMSArrayEvent.Delete();
968 fPadPedestalArrayEvent.Delete();
970 for ( Int_t i=0; i<72; ++i ){
972 fVTime0OffsetCounter[i]=0;
975 //_____________________________________________________________________
976 void AliTPCCalibPulser::ResetPad()
979 // Reset pad infos -- Should be called after a pad has been processed
981 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
988 //_____________________________________________________________________
989 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
992 // Merge reference histograms of sig to the current AliTPCCalibPulser
996 for (Int_t iSec=0; iSec<72; ++iSec){
997 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
998 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
999 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1003 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1004 TH2S *hRefQ = GetHistoQ(iSec);
1005 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1007 TH2S *hist = new TH2S(*hRefQmerge);
1008 hist->SetDirectory(0);
1009 fHistoQArray.AddAt(hist, iSec);
1011 hRefQmerge->SetDirectory(dir);
1014 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1015 TH2S *hRefT0 = GetHistoT0(iSec);
1016 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1018 TH2S *hist = new TH2S(*hRefT0merge);
1019 hist->SetDirectory(0);
1020 fHistoT0Array.AddAt(hist, iSec);
1022 hRefT0merge->SetDirectory(dir);
1024 if ( hRefRMSmerge ){
1025 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1026 TH2S *hRefRMS = GetHistoRMS(iSec);
1027 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1029 TH2S *hist = new TH2S(*hRefRMSmerge);
1030 hist->SetDirectory(0);
1031 fHistoRMSArray.AddAt(hist, iSec);
1033 hRefRMSmerge->SetDirectory(dir);
1038 //_____________________________________________________________________
1039 void AliTPCCalibPulser::Analyse()
1042 // Calculate calibration constants
1046 TVectorD paramT0(3);
1047 TVectorD paramRMS(3);
1048 TMatrixD dummy(3,3);
1050 for (Int_t iSec=0; iSec<72; ++iSec){
1051 TH2S *hT0 = GetHistoT0(iSec);
1052 if (!hT0 ) continue;
1054 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1055 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1056 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1057 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1059 TH2S *hQ = GetHistoQ(iSec);
1060 TH2S *hRMS = GetHistoRMS(iSec);
1062 Short_t *arrayhQ = hQ->GetArray();
1063 Short_t *arrayhT0 = hT0->GetArray();
1064 Short_t *arrayhRMS = hRMS->GetArray();
1066 UInt_t nChannels = fROC->GetNChannels(iSec);
1074 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1076 Float_t cogTime0 = -1000;
1077 Float_t cogQ = -1000;
1078 Float_t cogRMS = -1000;
1081 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1082 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1083 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1085 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,¶mQ,&dummy);
1086 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,¶mT0,&dummy);
1087 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,¶mRMS,&dummy);
1089 cogTime0 = paramT0[1];
1090 cogRMS = paramRMS[1];
1092 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1093 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1094 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1097 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1104 rocQ->SetValue(iChannel, cogQ*cogQ);
1105 rocT0->SetValue(iChannel, cogTime0);
1106 rocRMS->SetValue(iChannel, cogRMS);
1107 rocOut->SetValue(iChannel, cogOut);
1111 if ( fDebugLevel > 0 ){
1112 if ( !fDebugStreamer ) {
1114 TDirectory *backup = gDirectory;
1115 fDebugStreamer = new TTreeSRedirector("deb2.root");
1116 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1119 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1120 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1121 padc = pad-(fROC->GetNPads(iSec,row)/2);
1123 (*fDebugStreamer) << "DataEnd" <<
1124 "Sector=" << iSec <<
1128 "PadSec=" << iChannel <<
1130 "T0=" << cogTime0 <<
1137 delete fDebugStreamer;
1138 fDebugStreamer = 0x0;
1140 //_____________________________________________________________________
1141 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1144 // Write class to file
1153 option = "recreate";
1155 TDirectory *backup = gDirectory;
1156 TFile f(filename,option.Data());
1158 if ( !sDir.IsNull() ){
1159 f.mkdir(sDir.Data());
1165 if ( backup ) backup->cd();
1167 //_____________________________________________________________________
1168 //_________________________ Test Functions ___________________________
1169 //_____________________________________________________________________
1170 TObjArray* AliTPCCalibPulser::TestBinning()
1173 // Function to test the binning of the reference histograms
1174 // type: T0, Q or RMS
1175 // mode: 0 - number of filled bins per channel
1176 // 1 - number of empty bins between filled bins in one ROC
1177 // returns TObjArray with the test histograms type*2+mode:
1178 // position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1181 TObjArray *histArray = new TObjArray(6);
1182 const Char_t *type[] = {"T0","Q","RMS"};
1183 Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1185 for (Int_t itype = 0; itype<3; ++itype){
1186 for (Int_t imode=0; imode<2; ++imode){
1187 Int_t icount = itype*2+imode;
1188 histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1189 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1198 for (Int_t itype = 0; itype<3; ++itype){
1199 for (Int_t iSec=0; iSec<72; ++iSec){
1200 if ( itype == 0 ) hRef = GetHistoT0(iSec);
1201 if ( itype == 1 ) hRef = GetHistoQ(iSec);
1202 if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1203 if ( hRef == 0x0 ) continue;
1204 array = (hRef->GetArray());
1205 UInt_t nChannels = fROC->GetNChannels(iSec);
1208 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1210 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1213 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1214 if ( array[offset+iBin]>0 ) {
1216 if ( c1 && c2 ) nempty++;
1219 else if ( c1 ) c2 = 1;
1222 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1224 ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);