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() :
217 fIsZeroSuppressed(kFALSE),
219 fROC(AliTPCROC::Instance()),
221 fParam(new AliTPCParam),
230 fCalRocArrayOutliers(72),
234 fHMeanTimeSector(0x0),
235 fVMeanTimeSector(72),
236 fPadTimesArrayEvent(72),
238 fPadRMSArrayEvent(72),
239 fPadPedestalArrayEvent(72),
250 fVTime0OffsetCounter(72),
256 // AliTPCSignal default constructor
261 //_____________________________________________________________________
262 AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
264 fFirstTimeBin(sig.fFirstTimeBin),
265 fLastTimeBin(sig.fLastTimeBin),
266 fNbinsT0(sig.fNbinsT0),
267 fXminT0(sig.fXminT0),
268 fXmaxT0(sig.fXmaxT0),
269 fNbinsQ(sig.fNbinsQ),
272 fNbinsRMS(sig.fNbinsRMS),
273 fXminRMS(sig.fXminRMS),
274 fXmaxRMS(sig.fXmaxRMS),
275 fIsZeroSuppressed(sig.fIsZeroSuppressed),
277 fROC(AliTPCROC::Instance()),
279 fParam(new AliTPCParam),
288 fCalRocArrayOutliers(72),
292 fHMeanTimeSector(0x0),
293 fVMeanTimeSector(72),
294 fPadTimesArrayEvent(72),
296 fPadRMSArrayEvent(72),
297 fPadPedestalArrayEvent(72),
308 fVTime0OffsetCounter(72),
311 fDebugLevel(sig.fDebugLevel)
314 // AliTPCSignal default constructor
317 for (Int_t iSec = 0; iSec < 72; ++iSec){
318 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
319 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
320 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
321 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
323 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
324 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
325 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
327 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
328 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
329 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
330 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
333 TH2S *hNew = new TH2S(*hQ);
334 hNew->SetDirectory(0);
335 fHistoQArray.AddAt(hNew,iSec);
338 TH2S *hNew = new TH2S(*hT0);
339 hNew->SetDirectory(0);
340 fHistoQArray.AddAt(hNew,iSec);
343 TH2S *hNew = new TH2S(*hRMS);
344 hNew->SetDirectory(0);
345 fHistoQArray.AddAt(hNew,iSec);
347 fVMeanTimeSector[iSec]=sig.fVMeanTimeSector[iSec];
350 if ( sig.fHMeanTimeSector ) fHMeanTimeSector=(TH2F*)sig.fHMeanTimeSector->Clone();
353 AliTPCCalibPulser::AliTPCCalibPulser(const TMap *config) :
366 fIsZeroSuppressed(kFALSE),
368 fROC(AliTPCROC::Instance()),
370 fParam(new AliTPCParam),
379 fCalRocArrayOutliers(72),
383 fHMeanTimeSector(0x0),
384 fVMeanTimeSector(72),
385 fPadTimesArrayEvent(72),
387 fPadRMSArrayEvent(72),
388 fPadPedestalArrayEvent(72),
399 fVTime0OffsetCounter(72),
405 // This constructor uses a TMap for setting some parametes
407 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
408 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
409 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
410 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
411 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
412 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
413 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
414 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
415 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
416 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
417 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
418 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
422 //_____________________________________________________________________
423 AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
426 // assignment operator
428 if (&source == this) return *this;
429 new (this) AliTPCCalibPulser(source);
433 //_____________________________________________________________________
434 AliTPCCalibPulser::~AliTPCCalibPulser()
442 if ( fDebugStreamer) delete fDebugStreamer;
446 void AliTPCCalibPulser::Reset()
449 // Delete all information: Arrays, Histograms, CalRoc objects
451 fCalRocArrayT0.Delete();
452 fCalRocArrayQ.Delete();
453 fCalRocArrayRMS.Delete();
454 fCalRocArrayOutliers.Delete();
456 fHistoQArray.Delete();
457 fHistoT0Array.Delete();
458 fHistoRMSArray.Delete();
460 fPadTimesArrayEvent.Delete();
461 fPadQArrayEvent.Delete();
462 fPadRMSArrayEvent.Delete();
463 fPadPedestalArrayEvent.Delete();
465 if (fHMeanTimeSector) delete fHMeanTimeSector;
467 //_____________________________________________________________________
468 Int_t AliTPCCalibPulser::Update(const Int_t icsector,
471 const Int_t icTimeBin,
472 const Float_t csignal)
475 // Signal filling methode on the fly pedestal and time offset correction if necessary.
476 // no extra analysis necessary. Assumes knowledge of the signal shape!
477 // assumes that it is looped over consecutive time bins of one pad
480 if (icRow<0) return 0;
481 if (icPad<0) return 0;
482 if (icTimeBin<0) return 0;
483 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
485 if ( icRow<0 || icPad<0 ){
486 AliWarning("Wrong Pad or Row number, skipping!");
490 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
492 //init first pad and sector in this event
493 if ( fCurrentChannel == -1 ) {
494 fCurrentChannel = iChannel;
495 fCurrentSector = icsector;
500 //process last pad if we change to a new one
501 if ( iChannel != fCurrentChannel ){
503 fLastSector=fCurrentSector;
504 fCurrentChannel = iChannel;
505 fCurrentSector = icsector;
510 //fill signals for current pad
511 fPadSignal[icTimeBin]=csignal;
512 if ( csignal > fMaxPadSignal ){
513 fMaxPadSignal = csignal;
514 fMaxTimeBin = icTimeBin;
518 //_____________________________________________________________________
519 void AliTPCCalibPulser::FindPedestal(Float_t part)
522 // find pedestal and noise for the current pad. Use either database or
523 // truncated mean with part*100%
525 Bool_t noPedestal = kTRUE;;
526 if (fPedestalTPC&&fPadNoiseTPC){
527 //use pedestal database
528 //only load new pedestals if the sector has changed
529 if ( fCurrentSector!=fLastSector ){
530 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
531 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
534 if ( fPedestalROC&&fPadNoiseROC ){
535 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
536 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
542 //if we are not running with pedestal database, or for the current sector there is no information
543 //available, calculate the pedestal and noise on the fly
545 const Int_t kPedMax = 100; //maximum pedestal value
554 UShort_t histo[kPedMax];
555 memset(histo,0,kPedMax*sizeof(UShort_t));
557 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
558 padSignal = fPadSignal.GetMatrixArray()[i];
559 if (padSignal<=0) continue;
560 if (padSignal>max && i>10) {
564 if (padSignal>kPedMax-1) continue;
565 histo[Int_t(padSignal+0.5)]++;
569 for (Int_t i=1; i<kPedMax; ++i){
570 if (count1<count0*0.5) median=i;
575 // what if by chance histo[median] == 0 ?!?
576 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
578 for (Int_t idelta=1; idelta<10; ++idelta){
579 if (median-idelta<=0) continue;
580 if (median+idelta>kPedMax) continue;
581 if (count<part*count1){
582 count+=histo[median-idelta];
583 mean +=histo[median-idelta]*(median-idelta);
584 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
585 count+=histo[median+idelta];
586 mean +=histo[median+idelta]*(median+idelta);
587 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
594 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
599 fPadPedestal*=(Float_t)(!fIsZeroSuppressed);
601 //_____________________________________________________________________
602 void AliTPCCalibPulser::FindPulserSignal(TVectorD ¶m, Float_t &qSum)
605 // Find position, signal width and height of the CE signal (last signal)
606 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
607 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
610 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
611 Int_t cemaxpos = fMaxTimeBin;
612 Float_t ceSumThreshold = 10.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal sum
613 Float_t ceMaxThreshold = 5.*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal max
614 const Int_t kCemin = 2; // range for the analysis of the ce signal +- channels from the peak
615 const Int_t kCemax = 7;
622 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
623 if ( ceQmax<ceMaxThreshold ) return;
624 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
625 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
626 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
627 ceTime+=signal*(i+0.5);
628 ceRMS +=signal*(i+0.5)*(i+0.5);
633 if (ceQsum>ceSumThreshold) {
635 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
636 //only fill the Time0Offset if pad was not marked as an outlier!
638 //skip edge pads for calculating the mean time
639 if ( !IsEdgePad(fCurrentSector,fCurrentRow,fCurrentPad) ){
640 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
641 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
642 GetHistoTSec()->Fill(ceTime,fCurrentSector);
645 if ( !(fOutliers->GetCalROC(fCurrentSector)->GetValue(fCurrentChannel)) ){
646 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
647 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
651 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
652 // the pick-up signal should scale with the pad area. In addition
653 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
654 // ratio 2/3. The pad area we express in mm2 (factor 100). We normalise the signal
655 // to the OROC signal (factor 2/3 for the IROCs).
656 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow)*100;
657 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
671 //_____________________________________________________________________
672 void AliTPCCalibPulser::ProcessPad()
675 // Process data of current pad
681 FindPulserSignal(param, qSum);
683 Double_t meanT = param[1];
684 Double_t sigmaT = param[2];
687 //Fill Event T0 counter
688 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
691 // GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel ); //use linear scale, needs also a change in the Analyse funciton.
692 GetHistoQ(fCurrentSector,kTRUE)->Fill( qSum, fCurrentChannel );
695 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
698 //Fill debugging info
699 if ( fDebugLevel>0 ){
700 if ( fDebugLevel == 1 ){
701 if ( !fDebugStreamer ) {
703 TDirectory *backup = gDirectory;
704 fDebugStreamer = new TTreeSRedirector("debPulserPadSignals.root");
705 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
707 Int_t padc = fCurrentPad-(fROC->GetNPads(fCurrentSector,fCurrentRow)/2);
708 (*fDebugStreamer) << "PadSignals" <<
709 "Sector=" <<fCurrentSector<<
710 "Row=" <<fCurrentRow<<
711 "Pad=" <<fCurrentPad<<
713 "Channel="<<fCurrentChannel<<
716 "signal.=" <<&fPadSignal<<
719 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
720 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
721 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
726 //_____________________________________________________________________
727 void AliTPCCalibPulser::EndEvent()
730 // Process data of current event
733 //check if last pad has allready been processed, if not do so
734 if ( fMaxTimeBin>-1 ) ProcessPad();
736 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
737 for ( Int_t iSec = 0; iSec<72; ++iSec ){
738 TVectorF *vTimes = GetPadTimesEvent(iSec);
739 if ( !vTimes || fVTime0OffsetCounter[iSec]==0 ) continue;
740 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
741 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
742 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
744 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
745 //GetHistoT0(iSec,kTRUE)->Fill( time,iChannel );
749 if ( fDebugLevel>1 ){
750 if ( !fDebugStreamer ) {
752 TDirectory *backup = gDirectory;
753 fDebugStreamer = new TTreeSRedirector("deb2.root");
754 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
761 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
762 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
764 UInt_t channel=iChannel;
767 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
768 pad = channel-fROC->GetRowIndexes(sector)[row];
769 padc = pad-(fROC->GetNPads(sector,row)/2);
771 (*fDebugStreamer) << "DataPad" <<
772 // "Event=" << fEvent <<
773 "Sector="<< sector <<
777 "PadSec="<< channel <<
788 //_____________________________________________________________________
789 Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
792 // Event Processing loop - AliTPCRawStream
795 Bool_t withInput = kFALSE;
796 while ( rawStreamFast->NextDDL() ){
797 while ( rawStreamFast->NextChannel() ){
798 Int_t isector = rawStreamFast->GetSector(); // current sector
799 Int_t iRow = rawStreamFast->GetRow(); // current row
800 Int_t iPad = rawStreamFast->GetPad(); // current pad
802 while ( rawStreamFast->NextBunch() ){
803 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
804 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
805 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
806 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
807 Update(isector,iRow,iPad,iTimeBin+1,signal);
818 //_____________________________________________________________________
819 Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
822 // Event processing loop - AliRawReader
824 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
825 Bool_t res=ProcessEventFast(rawStreamFast);
826 delete rawStreamFast;
829 //_____________________________________________________________________
830 Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
833 // Event Processing loop - AliTPCRawStream
838 Bool_t withInput = kFALSE;
840 while (rawStream->Next()) {
841 Int_t isector = rawStream->GetSector(); // current sector
842 Int_t iRow = rawStream->GetRow(); // current row
843 Int_t iPad = rawStream->GetPad(); // current pad
844 Int_t iTimeBin = rawStream->GetTime(); // current time bin
845 Float_t signal = rawStream->GetSignal(); // current ADC signal
847 Update(isector,iRow,iPad,iTimeBin,signal);
855 //_____________________________________________________________________
856 Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
859 // Event processing loop - AliRawReader
863 AliTPCRawStream rawStream(rawReader, (AliAltroMapping**)fMapping);
865 rawReader->Select("TPC");
867 return ProcessEvent(&rawStream);
869 //_____________________________________________________________________
870 Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
873 // Event processing loop - date event
875 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
876 Bool_t result=ProcessEvent(rawReader);
881 //_____________________________________________________________________
882 TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
883 Int_t nbinsY, Float_t ymin, Float_t ymax,
884 const Char_t *type, Bool_t force)
887 // return pointer to Q histogram
888 // if force is true create a new histogram if it doesn't exist allready
890 if ( !force || arr->UncheckedAt(sector) )
891 return (TH2S*)arr->UncheckedAt(sector);
893 // if we are forced and histogram doesn't yes exist create it
894 Char_t name[255], title[255];
896 sprintf(name,"hCalib%s%.2d",type,sector);
897 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
899 // new histogram with Q calib information. One value for each pad!
900 TH2S* hist = new TH2S(name,title,
902 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
903 hist->SetDirectory(0);
904 arr->AddAt(hist,sector);
907 //_____________________________________________________________________
908 TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
911 // return pointer to T0 histogram
912 // if force is true create a new histogram if it doesn't exist allready
914 TObjArray *arr = &fHistoT0Array;
915 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
917 //_____________________________________________________________________
918 TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
921 // return pointer to Q histogram
922 // if force is true create a new histogram if it doesn't exist allready
924 TObjArray *arr = &fHistoQArray;
925 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
927 //_____________________________________________________________________
928 TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
931 // return pointer to Q histogram
932 // if force is true create a new histogram if it doesn't exist allready
934 TObjArray *arr = &fHistoRMSArray;
935 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
937 //_____________________________________________________________________
938 TH2F* AliTPCCalibPulser::GetHistoTSec()
941 // return the pointer to the abs time distribution per sector
942 // create it if it does not exist
944 if ( !fHMeanTimeSector ) //!!!if you change the binning here, you should also change it in the Analyse Function!!
945 fHMeanTimeSector = new TH2F("fHMeanTimeSector","Abs mean time per sector",
946 20*(fLastTimeBin-fFirstTimeBin), fFirstTimeBin, fLastTimeBin,
948 return fHMeanTimeSector;
950 //_____________________________________________________________________
951 TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
954 // return pointer to Pad Info from 'arr' for the current event and sector
955 // if force is true create it if it doesn't exist allready
957 if ( !force || arr->UncheckedAt(sector) )
958 return (TVectorF*)arr->UncheckedAt(sector);
960 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
961 arr->AddAt(vect,sector);
964 //_____________________________________________________________________
965 TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
968 // return pointer to Pad Times Array for the current event and sector
969 // if force is true create it if it doesn't exist allready
971 TObjArray *arr = &fPadTimesArrayEvent;
972 return GetPadInfoEvent(sector,arr,force);
974 //_____________________________________________________________________
975 TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
978 // return pointer to Pad Q Array for the current event and sector
979 // if force is true create it if it doesn't exist allready
980 // for debugging purposes only
983 TObjArray *arr = &fPadQArrayEvent;
984 return GetPadInfoEvent(sector,arr,force);
986 //_____________________________________________________________________
987 TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
990 // return pointer to Pad RMS Array for the current event and sector
991 // if force is true create it if it doesn't exist allready
992 // for debugging purposes only
994 TObjArray *arr = &fPadRMSArrayEvent;
995 return GetPadInfoEvent(sector,arr,force);
997 //_____________________________________________________________________
998 TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
1001 // return pointer to Pad RMS Array for the current event and sector
1002 // if force is true create it if it doesn't exist allready
1003 // for debugging purposes only
1005 TObjArray *arr = &fPadPedestalArrayEvent;
1006 return GetPadInfoEvent(sector,arr,force);
1008 //_____________________________________________________________________
1009 AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1012 // return pointer to ROC Calibration
1013 // if force is true create a new histogram if it doesn't exist allready
1015 if ( !force || arr->UncheckedAt(sector) )
1016 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1018 // if we are forced and histogram doesn't yes exist create it
1020 // new AliTPCCalROC for T0 information. One value for each pad!
1021 AliTPCCalROC *croc = new AliTPCCalROC(sector);
1022 arr->AddAt(croc,sector);
1025 //_____________________________________________________________________
1026 AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
1029 // return pointer to Carge ROC Calibration
1030 // if force is true create a new histogram if it doesn't exist allready
1032 TObjArray *arr = &fCalRocArrayT0;
1033 return GetCalRoc(sector, arr, force);
1035 //_____________________________________________________________________
1036 AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
1039 // return pointer to T0 ROC Calibration
1040 // if force is true create a new histogram if it doesn't exist allready
1042 TObjArray *arr = &fCalRocArrayQ;
1043 return GetCalRoc(sector, arr, force);
1045 //_____________________________________________________________________
1046 AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
1049 // return pointer to signal width ROC Calibration
1050 // if force is true create a new histogram if it doesn't exist allready
1052 TObjArray *arr = &fCalRocArrayRMS;
1053 return GetCalRoc(sector, arr, force);
1055 //_____________________________________________________________________
1056 AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
1059 // return pointer to Outliers
1060 // if force is true create a new histogram if it doesn't exist allready
1062 TObjArray *arr = &fCalRocArrayOutliers;
1063 return GetCalRoc(sector, arr, force);
1065 //_____________________________________________________________________
1066 void AliTPCCalibPulser::ResetEvent()
1069 // Reset global counters -- Should be called before each event is processed
1079 fPadTimesArrayEvent.Delete();
1081 fPadQArrayEvent.Delete();
1082 fPadRMSArrayEvent.Delete();
1083 fPadPedestalArrayEvent.Delete();
1085 for ( Int_t i=0; i<72; ++i ){
1087 fVTime0OffsetCounter[i]=0;
1090 //_____________________________________________________________________
1091 void AliTPCCalibPulser::ResetPad()
1094 // Reset pad infos -- Should be called after a pad has been processed
1096 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1103 //_____________________________________________________________________
1104 Bool_t AliTPCCalibPulser::IsEdgePad(Int_t sector, Int_t row, Int_t pad)
1107 // return true if pad is on the edge of a row
1110 Int_t edge2 = fROC->GetNPads(sector,row)-1;
1111 if ( pad == edge1 || pad == edge2 ) return kTRUE;
1115 //_____________________________________________________________________
1116 void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
1119 // Merge reference histograms of sig to the current AliTPCCalibPulser
1123 for (Int_t iSec=0; iSec<72; ++iSec){
1124 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
1125 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
1126 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
1130 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1131 TH2S *hRefQ = GetHistoQ(iSec);
1132 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1134 TH2S *hist = new TH2S(*hRefQmerge);
1135 hist->SetDirectory(0);
1136 fHistoQArray.AddAt(hist, iSec);
1138 hRefQmerge->SetDirectory(dir);
1141 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1142 TH2S *hRefT0 = GetHistoT0(iSec);
1143 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1145 TH2S *hist = new TH2S(*hRefT0merge);
1146 hist->SetDirectory(0);
1147 fHistoT0Array.AddAt(hist, iSec);
1149 hRefT0merge->SetDirectory(dir);
1151 if ( hRefRMSmerge ){
1152 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1153 TH2S *hRefRMS = GetHistoRMS(iSec);
1154 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1156 TH2S *hist = new TH2S(*hRefRMSmerge);
1157 hist->SetDirectory(0);
1158 fHistoRMSArray.AddAt(hist, iSec);
1160 hRefRMSmerge->SetDirectory(dir);
1164 if ( sig->fHMeanTimeSector ){
1165 TDirectory *dir = sig->fHMeanTimeSector->GetDirectory(); sig->fHMeanTimeSector->SetDirectory(0);
1166 if ( fHMeanTimeSector ) fHMeanTimeSector->Add(sig->fHMeanTimeSector);
1168 fHMeanTimeSector = new TH2F(*sig->fHMeanTimeSector);
1169 fHMeanTimeSector->SetDirectory(0);
1171 sig->fHMeanTimeSector->SetDirectory(dir);
1174 //_____________________________________________________________________
1175 void AliTPCCalibPulser::Analyse()
1178 // Calculate calibration constants
1182 TVectorD paramT0(3);
1183 TVectorD paramRMS(3);
1184 TMatrixD dummy(3,3);
1185 //calculate mean time for each sector and mean time for each side
1186 TH1F hMeanTsec("hMeanTsec","hMeanTsec",20*(fLastTimeBin-fFirstTimeBin),fFirstTimeBin,fLastTimeBin);
1187 fVMeanTimeSector.Zero();
1189 for (Int_t iSec=0; iSec<72; ++iSec){
1190 TH2S *hT0 = GetHistoT0(iSec);
1191 if (!hT0 ) continue;
1192 //calculate sector mean T
1193 if ( fHMeanTimeSector ){
1194 Int_t nbinsT = fHMeanTimeSector->GetNbinsX();
1195 Int_t offset = (nbinsT+2)*(iSec+1);
1196 Float_t *arrP=fHMeanTimeSector->GetArray()+offset;
1198 for ( Int_t i=0; i<nbinsT; i++ ) entries+=(Int_t)arrP[i+1];
1199 hMeanTsec.Set(nbinsT+2,arrP);
1200 hMeanTsec.SetEntries(entries);
1202 // truncated mean: remove lower 5% and upper 5%
1203 if ( entries>0 ) AliMathBase::TruncatedMean(&hMeanTsec,¶mT0,0.05,.95);
1204 fVMeanTimeSector[iSec]=paramT0[1];
1207 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1208 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1209 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1210 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1212 TH2S *hQ = GetHistoQ(iSec);
1213 TH2S *hRMS = GetHistoRMS(iSec);
1215 Short_t *arrayhQ = hQ->GetArray();
1216 Short_t *arrayhT0 = hT0->GetArray();
1217 Short_t *arrayhRMS = hRMS->GetArray();
1219 UInt_t nChannels = fROC->GetNChannels(iSec);
1220 Float_t meanTsec = fVMeanTimeSector[iSec];
1228 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1230 Float_t cogTime0 = -1000;
1231 Float_t cogQ = -1000;
1232 Float_t cogRMS = -1000;
1235 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1236 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1237 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1239 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,¶mQ,&dummy);
1240 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,¶mT0,&dummy);
1241 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,¶mRMS,&dummy);
1243 cogTime0 = paramT0[1];
1244 cogRMS = paramRMS[1];
1246 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1247 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1248 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1251 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1258 // rocQ->SetValue(iChannel, cogQ*cogQ); // changed to linear scale again
1259 rocQ->SetValue(iChannel, cogQ);
1260 rocT0->SetValue(iChannel, cogTime0+meanTsec); //offset by mean time of the sector
1261 rocRMS->SetValue(iChannel, cogRMS);
1262 rocOut->SetValue(iChannel, cogOut);
1266 if ( fDebugLevel > 2 ){
1267 if ( !fDebugStreamer ) {
1269 TDirectory *backup = gDirectory;
1270 fDebugStreamer = new TTreeSRedirector("deb2.root");
1271 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1274 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1275 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1276 padc = pad-(fROC->GetNPads(iSec,row)/2);
1278 (*fDebugStreamer) << "DataEnd" <<
1279 "Sector=" << iSec <<
1283 "PadSec=" << iChannel <<
1285 "T0=" << cogTime0 <<
1294 delete fDebugStreamer;
1295 fDebugStreamer = 0x0;
1297 //_____________________________________________________________________
1298 void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1301 // Write class to file
1310 option = "recreate";
1312 TDirectory *backup = gDirectory;
1313 TFile f(filename,option.Data());
1315 if ( !sDir.IsNull() ){
1316 f.mkdir(sDir.Data());
1322 if ( backup ) backup->cd();
1324 //_____________________________________________________________________
1325 //_________________________ Test Functions ___________________________
1326 //_____________________________________________________________________
1327 TObjArray* AliTPCCalibPulser::TestBinning()
1330 // Function to test the binning of the reference histograms
1331 // type: T0, Q or RMS
1332 // mode: 0 - number of filled bins per channel
1333 // 1 - number of empty bins between filled bins in one ROC
1334 // returns TObjArray with the test histograms type*2+mode:
1335 // position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1338 TObjArray *histArray = new TObjArray(6);
1339 const Char_t *type[] = {"T0","Q","RMS"};
1340 Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1342 for (Int_t itype = 0; itype<3; ++itype){
1343 for (Int_t imode=0; imode<2; ++imode){
1344 Int_t icount = itype*2+imode;
1345 histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1346 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1355 for (Int_t itype = 0; itype<3; ++itype){
1356 for (Int_t iSec=0; iSec<72; ++iSec){
1357 if ( itype == 0 ) hRef = GetHistoT0(iSec);
1358 if ( itype == 1 ) hRef = GetHistoQ(iSec);
1359 if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1360 if ( hRef == 0x0 ) continue;
1361 array = (hRef->GetArray());
1362 UInt_t nChannels = fROC->GetNChannels(iSec);
1365 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1367 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1370 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1371 if ( array[offset+iBin]>0 ) {
1373 if ( c1 && c2 ) nempty++;
1376 else if ( c1 ) c2 = 1;
1379 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1381 ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);