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>
172 #include <TVectorF.h>
175 #include <TDirectory.h>
181 #include "AliRawReader.h"
182 #include "AliRawReaderRoot.h"
183 #include "AliRawReaderDate.h"
184 #include "AliTPCRawStream.h"
185 #include "AliTPCRawStreamFast.h"
186 #include "AliTPCCalROC.h"
187 #include "AliTPCCalPad.h"
188 #include "AliTPCROC.h"
189 #include "AliTPCParam.h"
190 #include "AliTPCCalibPulser.h"
191 #include "AliTPCcalibDB.h"
192 #include "AliMathBase.h"
194 #include "TTreeStream.h"
202 ClassImp(AliTPCCalibPulser)
204 AliTPCCalibPulser::AliTPCCalibPulser() :
205 AliTPCCalibRawBase(),
217 fIsZeroSuppressed(kFALSE),
219 fParam(new AliTPCParam),
228 fCalRocArrayOutliers(72),
232 fHMeanTimeSector(0x0),
233 fVMeanTimeSector(72),
234 fPadTimesArrayEvent(72),
236 fPadRMSArrayEvent(72),
237 fPadPedestalArrayEvent(72),
248 fVTime0OffsetCounter(72)
251 // AliTPCSignal default constructor
253 SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
258 //_____________________________________________________________________
259 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
260 AliTPCCalibRawBase(sig),
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),
270 fPeakIntMinus(sig.fPeakIntMinus),
271 fPeakIntPlus(sig.fPeakIntPlus),
272 fIsZeroSuppressed(sig.fIsZeroSuppressed),
274 fParam(new AliTPCParam),
283 fCalRocArrayOutliers(72),
287 fHMeanTimeSector(0x0),
288 fVMeanTimeSector(72),
289 fPadTimesArrayEvent(72),
291 fPadRMSArrayEvent(72),
292 fPadPedestalArrayEvent(72),
303 fVTime0OffsetCounter(72)
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);
339 fVMeanTimeSector[iSec]=sig.fVMeanTimeSector[iSec];
342 if ( sig.fHMeanTimeSector ) fHMeanTimeSector=(TH2F*)sig.fHMeanTimeSector->Clone();
345 AliTPCCalibPulser::AliTPCCalibPulser(const TMap *config) :
346 AliTPCCalibRawBase(),
358 fIsZeroSuppressed(kFALSE),
360 fParam(new AliTPCParam),
369 fCalRocArrayOutliers(72),
373 fHMeanTimeSector(0x0),
374 fVMeanTimeSector(72),
375 fPadTimesArrayEvent(72),
377 fPadRMSArrayEvent(72),
378 fPadPedestalArrayEvent(72),
389 fVTime0OffsetCounter(72)
392 // This constructor uses a TMap for setting some parametes
394 SetNameTitle("AliTPCCalibPulser","AliTPCCalibPulser");
397 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
398 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
399 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
400 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
401 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
402 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
403 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
404 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
405 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
406 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
407 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
408 if (config->GetValue("PeakIntMinus")) fPeakIntMinus = (Int_t)((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atof();
409 if (config->GetValue("PeakIntPlus")) fPeakIntPlus = (Int_t)((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atof();
410 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
414 //_____________________________________________________________________
415 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
418 // assignment operator
420 if (&source == this) return *this;
421 new (this) AliTPCCalibPulser(source);
425 //_____________________________________________________________________
426 AliTPCCalibPulser::~AliTPCCalibPulser()
435 void AliTPCCalibPulser::Reset()
438 // Delete all information: Arrays, Histograms, CalRoc objects
440 fCalRocArrayT0.Delete();
441 fCalRocArrayQ.Delete();
442 fCalRocArrayRMS.Delete();
443 fCalRocArrayOutliers.Delete();
445 fHistoQArray.Delete();
446 fHistoT0Array.Delete();
447 fHistoRMSArray.Delete();
449 fPadTimesArrayEvent.Delete();
450 fPadQArrayEvent.Delete();
451 fPadRMSArrayEvent.Delete();
452 fPadPedestalArrayEvent.Delete();
454 if (fHMeanTimeSector) delete fHMeanTimeSector;
456 //_____________________________________________________________________
457 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
460 const Int_t icTimeBin,
461 const Float_t csignal)
464 // Signal filling methode on the fly pedestal and time offset correction if necessary.
465 // no extra analysis necessary. Assumes knowledge of the signal shape!
466 // assumes that it is looped over consecutive time bins of one pad
469 if (icRow<0) return 0;
470 if (icPad<0) return 0;
471 if (icTimeBin<0) return 0;
472 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
474 if ( icRow<0 || icPad<0 ){
475 AliWarning("Wrong Pad or Row number, skipping!");
479 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
481 //init first pad and sector in this event
482 if ( fCurrentChannel == -1 ) {
483 fCurrentChannel = iChannel;
484 fCurrentSector = icsector;
489 //process last pad if we change to a new one
490 if ( iChannel != fCurrentChannel ){
492 fLastSector=fCurrentSector;
493 fCurrentChannel = iChannel;
494 fCurrentSector = icsector;
499 //fill signals for current pad
500 fPadSignal[icTimeBin]=csignal;
501 if ( csignal > fMaxPadSignal ){
502 fMaxPadSignal = csignal;
503 fMaxTimeBin = icTimeBin;
507 //_____________________________________________________________________
508 void AliTPCCalibPulser::FindPedestal(Float_t part)
511 // find pedestal and noise for the current pad. Use either database or
512 // truncated mean with part*100%
514 Bool_t noPedestal = kTRUE;;
515 if (fPedestalTPC&&fPadNoiseTPC){
516 //use pedestal database
517 //only load new pedestals if the sector has changed
518 if ( fCurrentSector!=fLastSector ){
519 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
520 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
523 if ( fPedestalROC&&fPadNoiseROC ){
524 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
525 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
531 //if we are not running with pedestal database, or for the current sector there is no information
532 //available, calculate the pedestal and noise on the fly
534 const Int_t kPedMax = 100; //maximum pedestal value
543 UShort_t histo[kPedMax];
544 memset(histo,0,kPedMax*sizeof(UShort_t));
546 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
547 padSignal = fPadSignal.GetMatrixArray()[i];
548 if (padSignal<=0) continue;
549 if (padSignal>max && i>10) {
553 if (padSignal>kPedMax-1) continue;
554 histo[Int_t(padSignal+0.5)]++;
558 for (Int_t i=1; i<kPedMax; ++i){
559 if (count1<count0*0.5) median=i;
564 // what if by chance histo[median] == 0 ?!?
565 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
567 for (Int_t idelta=1; idelta<10; ++idelta){
568 if (median-idelta<=0) continue;
569 if (median+idelta>kPedMax) continue;
570 if (count<part*count1){
571 count+=histo[median-idelta];
572 mean +=histo[median-idelta]*(median-idelta);
573 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
574 count+=histo[median+idelta];
575 mean +=histo[median+idelta]*(median+idelta);
576 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
583 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
588 fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
590 //_____________________________________________________________________
591 void AliTPCCalibPulser::FindPulserSignal(TVectorD ¶m, Float_t &qSum)
594 // Find position, signal width and height of the CE signal (last signal)
595 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
596 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
599 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
600 Int_t cemaxpos = fMaxTimeBin;
601 Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal sum
602 Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal max
603 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
604 const Int_t kCemax = fPeakIntPlus;
611 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
612 if ( ceQmax<ceMaxThreshold ) return;
613 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
614 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
615 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
616 ceTime+=signal*(i+0.5);
617 ceRMS +=signal*(i+0.5)*(i+0.5);
622 if (ceQsum>ceSumThreshold) {
624 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
625 ceTime-=GetL1PhaseTB();
626 //only fill the Time0Offset if pad was not marked as an outlier!
628 //skip edge pads for calculating the mean time
629 if ( !IsEdgePad(fCurrentSector,fCurrentRow,fCurrentPad) ){
630 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
631 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
632 GetHistoTSec()->Fill(ceTime,fCurrentSector);
635 if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
636 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
637 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
641 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
642 // the pick-up signal should scale with the pad area. In addition
643 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
644 // ratio 2/3. The pad area we express in mm2 (factor 100). We normalise the signal
645 // to the OROC signal (factor 2/3 for the IROCs).
646 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow)*100;
647 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
661 //_____________________________________________________________________
662 void AliTPCCalibPulser::ProcessPad()
665 // Process data of current pad
671 FindPulserSignal(param, qSum);
673 Double_t meanT = param[1];
674 Double_t sigmaT = param[2];
677 //Fill Event T0 counter
678 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
681 // GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel ); //use linear scale, needs also a change in the Analyse funciton.
682 GetHistoQ(fCurrentSector,kTRUE)->Fill( qSum, fCurrentChannel );
685 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
688 //Fill debugging info
689 if ( GetStreamLevel()>0 ){
690 TTreeSRedirector *streamer=GetDebugStreamer();
691 if ( GetStreamLevel() == 1 ){
693 Int_t padc = fCurrentPad-(fROC->GetNPads(fCurrentSector,fCurrentRow)/2);
694 (*streamer) << "PadSignals" <<
695 "Event=" <<fNevents <<
696 "Sector=" <<fCurrentSector<<
697 "Row=" <<fCurrentRow<<
698 "Pad=" <<fCurrentPad<<
700 "Channel="<<fCurrentChannel<<
703 "signal.=" <<&fPadSignal<<
707 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
708 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
709 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
714 //_____________________________________________________________________
715 void AliTPCCalibPulser::EndEvent()
718 // Process data of current event
721 //check if last pad has allready been processed, if not do so
722 if ( fMaxTimeBin>-1 ) ProcessPad();
724 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
725 for ( Int_t iSec = 0; iSec<72; ++iSec ){
726 TVectorF *vTimes = GetPadTimesEvent(iSec);
727 if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
728 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
729 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
730 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
732 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
733 //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
737 if ( GetStreamLevel()>1 ){
738 TTreeSRedirector *streamer=GetDebugStreamer();
744 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
745 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
747 UInt_t channel=iChannel;
750 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
751 pad = channel-fROC->GetRowIndexes(sector)[row];
752 padc = pad-(fROC->GetNPads(sector,row)/2);
754 (*streamer) << "DataPad" <<
755 "Event=" << fNevents <<
756 "Sector="<< sector <<
760 "PadSec="<< channel <<
773 //_____________________________________________________________________
774 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
775 Int_t nbinsY, Float_t ymin, Float_t ymax,
776 const Char_t *type, Bool_t force)
779 // return pointer to Q histogram
780 // if force is true create a new histogram if it doesn't exist allready
782 if ( !force || arr->UncheckedAt(sector) )
783 return (TH2S*)arr->UncheckedAt(sector);
785 // if we are forced and histogram doesn't yes exist create it
786 // new histogram with Q calib information. One value for each pad!
787 TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
789 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
790 hist->SetDirectory(0);
791 arr->AddAt(hist,sector);
794 //_____________________________________________________________________
795 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
798 // return pointer to T0 histogram
799 // if force is true create a new histogram if it doesn't exist allready
801 TObjArray *arr = &fHistoT0Array;
802 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
804 //_____________________________________________________________________
805 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
808 // return pointer to Q histogram
809 // if force is true create a new histogram if it doesn't exist allready
811 TObjArray *arr = &fHistoQArray;
812 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
814 //_____________________________________________________________________
815 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
818 // return pointer to Q histogram
819 // if force is true create a new histogram if it doesn't exist allready
821 TObjArray *arr = &fHistoRMSArray;
822 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
824 //_____________________________________________________________________
825 TH2F* AliTPCCalibPulser::GetHistoTSec()
828 // return the pointer to the abs time distribution per sector
829 // create it if it does not exist
831 if ( !fHMeanTimeSector ) //!!!if you change the binning here, you should also change it in the Analyse Function!!
832 fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
833 20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
835 return fHMeanTimeSector;
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
966 fPadTimesArrayEvent.Delete();
968 fPadQArrayEvent.Delete();
969 fPadRMSArrayEvent.Delete();
970 fPadPedestalArrayEvent.Delete();
972 for ( Int_t i=0; i<72; ++i ){
974 fVTime0OffsetCounter[i]=0;
977 //_____________________________________________________________________
978 void AliTPCCalibPulser::ResetPad()
981 // Reset pad infos -- Should be called after a pad has been processed
983 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
990 //_____________________________________________________________________
991 Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
994 // return true if pad is on the edge of a row
997 Int_t edge2 = fROC->GetNPads(sector,row)-1;
998 if ( pad == edge1 || pad == edge2 ) return kTRUE;
1002 //_____________________________________________________________________
1003 void AliTPCCalibPulser::Merge(AliTPCCalibPulser * const sig)
1006 // Merge reference histograms of sig to the current AliTPCCalibPulser
1011 for (Int_t iSec=0; iSec<72; ++iSec){
1012 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
1013 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
1014 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1018 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1019 TH2S *hRefQ = GetHistoQ(iSec);
1020 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1022 TH2S *hist = new TH2S(*hRefQmerge);
1023 hist->SetDirectory(0);
1024 fHistoQArray.AddAt(hist, iSec);
1026 hRefQmerge->SetDirectory(dir);
1029 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1030 TH2S *hRefT0 = GetHistoT0(iSec);
1031 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1033 TH2S *hist = new TH2S(*hRefT0merge);
1034 hist->SetDirectory(0);
1035 fHistoT0Array.AddAt(hist, iSec);
1037 hRefT0merge->SetDirectory(dir);
1039 if ( hRefRMSmerge ){
1040 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1041 TH2S *hRefRMS = GetHistoRMS(iSec);
1042 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1044 TH2S *hist = new TH2S(*hRefRMSmerge);
1045 hist->SetDirectory(0);
1046 fHistoRMSArray.AddAt(hist, iSec);
1048 hRefRMSmerge->SetDirectory(dir);
1052 if ( sig->fHMeanTimeSector ){
1053 TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1054 if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1056 fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1057 fHMeanTimeSector->SetDirectory(0);
1059 sig->fHMeanTimeSector->SetDirectory(dir);
1064 //_____________________________________________________________________
1065 Long64_t AliTPCCalibPulser::Merge(TCollection * const list)
1068 // Merge all objects of this type in list
1074 AliTPCCalibPulser *ce=0;
1077 while ( (o=next()) ){
1078 ce=dynamic_cast<AliTPCCalibPulser*>(o);
1088 //_____________________________________________________________________
1089 void AliTPCCalibPulser::Analyse()
1092 // Calculate calibration constants
1096 TVectorD paramT0(3);
1097 TVectorD paramRMS(3);
1098 TMatrixD dummy(3,3);
1099 //calculate mean time for each sector and mean time for each side
1100 TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1101 fVMeanTimeSector.Zero();
1103 for (Int_t iSec=0; iSec<72; ++iSec){
1104 TH2S *hT0 = GetHistoT0(iSec);
1105 if (!hT0 ) continue;
1106 //calculate sector mean T
1107 if ( fHMeanTimeSector ){
1108 Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1109 Int_t offset = (nbinsT+2)*(iSec+1);
1110 Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1112 for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1113 hMeanTsec.Set(nbinsT+2,arrP);
1114 hMeanTsec.SetEntries(entries);
1116 // truncated mean: remove lower 5% and upper 5%
1117 if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,¶mT0,0.05,.95);
1118 fVMeanTimeSector[iSec]=paramT0[1];
1121 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1122 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1123 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1124 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1126 TH2S *hQ = GetHistoQ(iSec);
1127 TH2S *hRMS = GetHistoRMS(iSec);
1129 Short_t *arrayhQ = hQ->GetArray();
1130 Short_t *arrayhT0 = hT0->GetArray();
1131 Short_t *arrayhRMS = hRMS->GetArray();
1133 UInt_t nChannels = fROC->GetNChannels(iSec);
1134 Float_t meanTsec = fVMeanTimeSector[iSec];
1142 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1144 Float_t cogTime0 = -1000;
1145 Float_t cogQ = -1000;
1146 Float_t cogRMS = -1000;
1149 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1150 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1151 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1153 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,¶mQ,&dummy);
1154 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,¶mT0,&dummy);
1155 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,¶mRMS,&dummy);
1157 cogTime0 = paramT0[1];
1158 cogRMS = paramRMS[1];
1160 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1161 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1162 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1165 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1172 // rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1173 rocQ->SetValue(iChannel, cogQ);
1174 rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1175 rocRMS->SetValue(iChannel, cogRMS);
1176 rocOut->SetValue(iChannel, cogOut);
1178 if ( GetStreamLevel() > 2 ){
1179 TTreeSRedirector *streamer=GetDebugStreamer();
1181 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1182 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1183 padc = pad-(fROC->GetNPads(iSec,row)/2);
1185 (*streamer) << "DataEnd" <<
1186 "Sector=" << iSec <<
1190 "PadSec=" << iChannel <<
1192 "T0=" << cogTime0 <<
1203 //_____________________________________________________________________
1204 //_________________________ Test Functions ___________________________
1205 //_____________________________________________________________________
1206 TObjArray* AliTPCCalibPulser::TestBinning()
1209 // Function to test the binning of the reference histograms
1210 // type: T0, Q or RMS
1211 // mode: 0 - number of filled bins per channel
1212 // 1 - number of empty bins between filled bins in one ROC
1213 // returns TObjArray with the test histograms type*2+mode:
1214 // position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1217 TObjArray *histArray = new TObjArray(6);
1218 const Char_t *type[] = {"T0","Q","RMS"};
1219 Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1221 for (Int_t itype = 0; itype<3; ++itype){
1222 for (Int_t imode=0; imode<2; ++imode){
1223 Int_t icount = itype*2+imode;
1224 histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1225 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1234 for (Int_t itype = 0; itype<3; ++itype){
1235 for (Int_t iSec=0; iSec<72; ++iSec){
1236 if ( itype == 0 ) hRef = GetHistoT0(iSec);
1237 if ( itype == 1 ) hRef = GetHistoQ(iSec);
1238 if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1239 if ( hRef == 0x0 ) continue;
1240 array = (hRef->GetArray());
1241 UInt_t nChannels = fROC->GetNChannels(iSec);
1244 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1246 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1249 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1250 if ( array[offset+iBin]>0 ) {
1252 if ( c1 && c2 ) nempty++;
1255 else if ( c1 ) c2 = 1;
1258 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1260 ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);