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 Central Electrode calibration //
22 // Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
24 ////////////////////////////////////////////////////////////////////////////////////////
27 // *************************************************************************************
28 // * Class Description *
29 // *************************************************************************************
32 <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
33 using laser runs.</h4>
35 The information retrieved is
36 <ul style="list-style-type: square;">
37 <li>Time arrival from the CE</li>
43 <ol style="list-style-type: upper-roman;">
44 <li><a href="#working">Working principle</a></li>
45 <li><a href="#user">User interface for filling data</a></li>
46 <li><a href="#info">Stored information</a></li>
49 <h3><a name="working">I. Working principle</a></h3>
51 <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
52 (see below). These in the end call the Update(...) function.</h4>
54 <ul style="list-style-type: square;">
55 <li>the Update(...) function:<br />
56 In this function the array fPadSignal is filled with the adc signals between the specified range
57 fFirstTimeBin and fLastTimeBin for the current pad.
58 before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
61 <ul style="list-style-type: square;">
62 <li>the ProcessPad() function:</li>
63 <ol style="list-style-type: decimal;">
64 <li>Find Pedestal and Noise information</li>
65 <ul style="list-style-type: square;">
66 <li>use database information which has to be set by calling<br />
67 SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
68 <li>if no information from the pedestal data base
69 is available the informaion is calculated on the fly
70 ( see FindPedestal() function )</li>
72 <li>Find local maxima of the pad signal</li>
73 <ul style="list-style-type: square;">
74 <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
75 have been observed ( see FindLocalMaxima(...) )</li>
77 <li>Find the CE signal information</li>
78 <ul style="list-style-type: square;">
79 <li>to find the position of the CE signal the Tmean information from the previos event is used
80 as the CE signal the local maximum closest to this Tmean is identified</li>
81 <li>calculate mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
82 the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
84 <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
85 <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
90 <h4>At the end of each event the EndEvent() function is called</h4>
92 <ul style="list-style-type: square;">
93 <li>the EndEvent() function:</li>
94 <ul style="list-style-type: square;">
95 <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
96 This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
97 <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
98 <li>calculate Mean Q</li>
99 <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
103 <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
104 <ul style="list-style-type: square;">
105 <li>the Analyse() function:</li>
106 <ul style="list-style-type: square;">
107 <li>calculate the mean values of T0, RMS, Q for each pad, using
108 the AliMathBase::GetCOG(...) function</li>
109 <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
110 (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
114 <h3><a name="user">II. User interface for filling data</a></h3>
116 <h4>To Fill information one of the following functions can be used:</h4>
118 <ul style="list-style-type: none;">
119 <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
120 <ul style="list-style-type: square;">
121 <li>process Date event</li>
122 <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
126 <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
127 <ul style="list-style-type: square;">
128 <li>process AliRawReader event</li>
129 <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
133 <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
134 <ul style="list-style-type: square;">
135 <li>process event from AliTPCRawStream</li>
136 <li>call Update function for signal filling</li>
140 <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
141 iPad, const Int_t iTimeBin, const Float_t signal);</li>
142 <ul style="list-style-type: square;">
143 <li>directly fill signal information (sector, row, pad, time bin, pad)
144 to the reference histograms</li>
148 <h4>It is also possible to merge two independently taken calibrations using the function</h4>
150 <ul style="list-style-type: none;">
151 <li> void Merge(AliTPCCalibSignal *sig)</li>
152 <ul style="list-style-type: square;">
153 <li>copy histograms in 'sig' if they do not exist in this instance</li>
154 <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
155 <li>After merging call Analyse again!</li>
160 <h4>example: filling data using root raw data:</h4>
162 void fillCE(Char_t *filename)
164 rawReader = new AliRawReaderRoot(fileName);
165 if ( !rawReader ) return;
166 AliTPCCalibCE *calib = new AliTPCCalibCE;
167 while (rawReader->NextEvent()){
168 calib->ProcessEvent(rawReader);
171 calib->DumpToFile("CEData.root");
177 <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
179 <h4><a name="info:stored">III.1 Stored information</a></h4>
180 <ul style="list-style-type: none;">
182 <ul style="list-style-type: none;">
183 <li>For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum])
184 is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
185 stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
189 <li>Calibration Data:</li>
190 <ul style="list-style-type: none;">
191 <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
192 the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
193 fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
194 is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
198 <li>For each event the following information is stored:</li>
200 <ul style="list-style-type: square;">
201 <li>event time ( TVectorD fVEventTime )</li>
202 <li>event id ( TVectorD fVEventNumber )</li>
204 <li>mean arrival time for each ROC ( TObjArray fTMeanArrayEvent )</li>
205 <li>mean Q for each ROC ( TObjArray fQMeanArrayEvent )</li>
206 <li>parameters of a plane fit for each ROC ( TObjArray fParamArrayEventPol1 )</li>
207 <li>parameters of a 2D parabola fit for each ROC ( TObjArray fParamArrayEventPol2 )</li>
211 <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
212 <ul style="list-style-type: none;">
213 <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
214 <ul style="list-style-type: square;">
215 <li>TH2F *GetHistoT0(Int_t sector);</li>
216 <li>TH2F *GetHistoRMS(Int_t sector);</li>
217 <li>TH2F *GetHistoQ(Int_t sector);</li>
221 <li>Accessing the calibration storage objects:</li>
222 <ul style="list-style-type: square;">
223 <li>AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values</li>
224 <li>AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values</li>
225 <li>AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values</li>
229 <li>Accessin the event by event information:</li>
230 <ul style="list-style-type: square;">
231 <li>The event by event information can be displayed using the</li>
232 <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
233 <li>which creates a graph from the specified variables</li>
237 <h4>example for visualisation:</h4>
239 //if the file "CEData.root" was created using the above example one could do the following:
240 TFile fileCE("CEData.root")
241 AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
242 ce->GetCalRocT0(0)->Draw("colz");
243 ce->GetCalRocRMS(0)->Draw("colz");
245 //or use the AliTPCCalPad functionality:
246 AliTPCCalPad padT0(ped->GetCalPadT0());
247 AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
248 padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
249 padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
251 //display event by event information:
252 //Draw mean arrival time as a function of the event time for oroc sector A00
253 ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
254 //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
255 ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");
258 //////////////////////////////////////////////////////////////////////////////////////
262 #include <TObjArray.h>
266 #include <TVectorF.h>
267 #include <TVectorD.h>
268 #include <TMatrixD.h>
273 #include <TDirectory.h>
278 #include "AliRawReader.h"
279 #include "AliRawReaderRoot.h"
280 #include "AliRawReaderDate.h"
281 #include "AliRawEventHeaderBase.h"
282 #include "AliTPCRawStream.h"
283 #include "AliTPCRawStreamFast.h"
284 #include "AliTPCcalibDB.h"
285 #include "AliTPCCalROC.h"
286 #include "AliTPCCalPad.h"
287 #include "AliTPCROC.h"
288 #include "AliTPCParam.h"
289 #include "AliTPCCalibCE.h"
290 #include "AliMathBase.h"
291 #include "TTreeStream.h"
295 ClassImp(AliTPCCalibCE)
298 AliTPCCalibCE::AliTPCCalibCE() :
312 fOldRCUformat(kTRUE),
313 fROC(AliTPCROC::Instance()),
314 fParam(new AliTPCParam),
322 fCalRocArrayOutliers(72),
327 fParamArrayEventPol1(72),
328 fParamArrayEventPol2(72),
329 fTMeanArrayEvent(72),
330 fQMeanArrayEvent(72),
338 fPadTimesArrayEvent(72),
340 fPadRMSArrayEvent(72),
341 fPadPedestalArrayEvent(72),
351 fVTime0OffsetCounter(72),
359 // AliTPCSignal default constructor
361 // fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
363 //_____________________________________________________________________
364 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
366 fFirstTimeBin(sig.fFirstTimeBin),
367 fLastTimeBin(sig.fLastTimeBin),
368 fNbinsT0(sig.fNbinsT0),
369 fXminT0(sig.fXminT0),
370 fXmaxT0(sig.fXmaxT0),
371 fNbinsQ(sig.fNbinsQ),
374 fNbinsRMS(sig.fNbinsRMS),
375 fXminRMS(sig.fXminRMS),
376 fXmaxRMS(sig.fXmaxRMS),
378 fOldRCUformat(kTRUE),
379 fROC(AliTPCROC::Instance()),
380 fParam(new AliTPCParam),
388 fCalRocArrayOutliers(72),
393 fParamArrayEventPol1(72),
394 fParamArrayEventPol2(72),
395 fTMeanArrayEvent(72),
396 fQMeanArrayEvent(72),
399 fNevents(sig.fNevents),
404 fPadTimesArrayEvent(72),
406 fPadRMSArrayEvent(72),
407 fPadPedestalArrayEvent(72),
417 fVTime0OffsetCounter(72),
422 fDebugLevel(sig.fDebugLevel)
425 // AliTPCSignal copy constructor
428 for (Int_t iSec = 0; iSec < 72; ++iSec){
429 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
430 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
431 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
432 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
434 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
435 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
436 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
438 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
439 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
440 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
441 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
444 // TDirectory *dir = hQ->GetDirectory();
445 // hQ->SetDirectory(0);
446 TH2S *hNew = new TH2S(*hQ);
447 hNew->SetDirectory(0);
448 fHistoQArray.AddAt(hNew,iSec);
449 // hQ->SetDirectory(dir);
452 // TDirectory *dir = hT0->GetDirectory();
453 // hT0->SetDirectory(0);
454 TH2S *hNew = new TH2S(*hT0);
455 hNew->SetDirectory(0);
456 fHistoT0Array.AddAt(hNew,iSec);
457 // hT0->SetDirectory(dir);
460 // TDirectory *dir = hRMS->GetDirectory();
461 // hRMS->SetDirectory(0);
462 TH2S *hNew = new TH2S(*hRMS);
463 hNew->SetDirectory(0);
464 fHistoRMSArray.AddAt(hNew,iSec);
465 // hRMS->SetDirectory(dir);
469 //copy fit parameters event by event
471 for (Int_t iSec=0; iSec<72; ++iSec){
472 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
474 TObjArray *arrEvents = new TObjArray(arr->GetSize());
475 fParamArrayEventPol1.AddAt(arrEvents, iSec);
476 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
477 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
478 arrEvents->AddAt(new TVectorD(*vec),iEvent);
481 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
483 TObjArray *arrEvents = new TObjArray(arr->GetSize());
484 fParamArrayEventPol2.AddAt(arrEvents, iSec);
485 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
486 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
487 arrEvents->AddAt(new TVectorD(*vec),iEvent);
490 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
491 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
493 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
495 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
499 fVEventTime.ResizeTo(sig.fVEventTime);
500 fVEventNumber.ResizeTo(sig.fVEventNumber);
501 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
502 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
505 //_____________________________________________________________________
506 AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
509 // assignment operator
511 if (&source == this) return *this;
512 new (this) AliTPCCalibCE(source);
516 //_____________________________________________________________________
517 AliTPCCalibCE::~AliTPCCalibCE()
523 fCalRocArrayT0.Delete();
524 fCalRocArrayQ.Delete();
525 fCalRocArrayRMS.Delete();
526 fCalRocArrayOutliers.Delete();
528 fHistoQArray.Delete();
529 fHistoT0Array.Delete();
530 fHistoRMSArray.Delete();
532 fHistoTmean.Delete();
534 fParamArrayEventPol1.Delete();
535 fParamArrayEventPol2.Delete();
536 fTMeanArrayEvent.Delete();
537 fQMeanArrayEvent.Delete();
539 fPadTimesArrayEvent.Delete();
540 fPadQArrayEvent.Delete();
541 fPadRMSArrayEvent.Delete();
542 fPadPedestalArrayEvent.Delete();
544 if ( fDebugStreamer) delete fDebugStreamer;
545 // if ( fHTime0 ) delete fHTime0;
549 //_____________________________________________________________________
550 Int_t AliTPCCalibCE::Update(const Int_t icsector,
553 const Int_t icTimeBin,
554 const Float_t csignal)
557 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
558 // no extra analysis necessary. Assumes knowledge of the signal shape!
559 // assumes that it is looped over consecutive time bins of one pad
561 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
563 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
565 //init first pad and sector in this event
566 if ( fCurrentChannel == -1 ) {
567 fCurrentChannel = iChannel;
568 fCurrentSector = icsector;
572 //process last pad if we change to a new one
573 if ( iChannel != fCurrentChannel ){
575 fCurrentChannel = iChannel;
576 fCurrentSector = icsector;
580 //fill signals for current pad
581 fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
582 if ( csignal > fMaxPadSignal ){
583 fMaxPadSignal = csignal;
584 fMaxTimeBin = icTimeBin;
588 //_____________________________________________________________________
589 void AliTPCCalibCE::FindPedestal(Float_t part)
592 // find pedestal and noise for the current pad. Use either database or
593 // truncated mean with part*100%
595 Bool_t noPedestal = kTRUE;
597 //use pedestal database if set
598 if (fPedestalTPC&&fPadNoiseTPC){
599 //only load new pedestals if the sector has changed
600 if ( fCurrentSector!=fLastSector ){
601 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
602 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
603 fLastSector=fCurrentSector;
606 if ( fPedestalROC&&fPadNoiseROC ){
607 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
608 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
614 //if we are not running with pedestal database, or for the current sector there is no information
615 //available, calculate the pedestal and noise on the fly
617 const Int_t kPedMax = 100; //maximum pedestal value
626 UShort_t histo[kPedMax];
627 memset(histo,0,kPedMax*sizeof(UShort_t));
629 //fill pedestal histogram
630 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
631 padSignal = fPadSignal.GetMatrixArray()[i];
632 if (padSignal<=0) continue;
633 if (padSignal>max && i>10) {
637 if (padSignal>kPedMax-1) continue;
638 histo[int(padSignal+0.5)]++;
642 for (Int_t i=1; i<kPedMax; ++i){
643 if (count1<count0*0.5) median=i;
648 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
650 for (Int_t idelta=1; idelta<10; ++idelta){
651 if (median-idelta<=0) continue;
652 if (median+idelta>kPedMax) continue;
653 if (count<part*count1){
654 count+=histo[median-idelta];
655 mean +=histo[median-idelta]*(median-idelta);
656 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
657 count+=histo[median+idelta];
658 mean +=histo[median+idelta]*(median+idelta);
659 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
666 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
672 //_____________________________________________________________________
673 void AliTPCCalibCE::FindCESignal(TVectorD ¶m, Float_t &qSum, const TVectorF maxima)
676 // Find position, signal width and height of the CE signal (last signal)
677 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
678 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
681 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
683 Float_t ceSumThreshold = 8.*fPadNoise; // threshold for the signal sum
684 const Int_t kCemin = 4; // range for the analysis of the ce signal +- channels from the peak
685 const Int_t kCemax = 7;
687 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
689 // find maximum closest to the sector mean from the last event
690 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
691 // get sector mean of last event
692 Float_t tmean = (*GetTMeanEvents(fCurrentSector))[fNevents-1];
693 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
694 minDist = tmean-maxima[imax];
695 cemaxpos = (Int_t)maxima[imax];
700 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
701 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
702 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
703 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
705 ceTime+=signal*(i+0.5);
706 ceRMS +=signal*(i+0.5)*(i+0.5);
712 if (ceQmax&&ceQsum>ceSumThreshold) {
714 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
715 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
716 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
718 //Normalise Q to pad area of irocs
719 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
722 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
723 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
735 //_____________________________________________________________________
736 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
739 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
740 // and 'tplus' timebins after 'pos'
742 if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
743 for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
744 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
745 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
746 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
748 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
752 //_____________________________________________________________________
753 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
756 // Find local maxima on the pad signal and Histogram them
758 Float_t ceThreshold = 5.*fPadNoise; // threshold for the signal
762 for (Int_t i=fLastTimeBin-tplus-1; i>=fFirstTimeBin+tminus; --i){
763 if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,tminus,tplus) ){
764 if (count<maxima.GetNrows()){
765 maxima.GetMatrixArray()[count++]=i;
766 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
771 //_____________________________________________________________________
772 void AliTPCCalibCE::ProcessPad()
775 // Process data of current pad
779 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
780 // + central electrode and possibly post peaks from the CE signal
781 // however if we are on a high noise pad a lot more peaks due to the noise might occur
782 FindLocalMaxima(maxima);
783 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
789 FindCESignal(param, qSum, maxima);
791 Double_t meanT = param[1];
792 Double_t sigmaT = param[2];
794 //Fill Event T0 counter
795 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
798 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
801 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
804 //Fill debugging info
805 if ( fDebugLevel>0 ){
806 (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
807 (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
808 (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
813 //_____________________________________________________________________
814 void AliTPCCalibCE::EndEvent()
817 // Process data of current pad
818 // The Functions 'SetTimeStamp' and 'SetRunNumber' should be called
819 // before the EndEvent function to set the event timestamp and number!!!
820 // This is automatically done if the ProcessEvent(AliRawReader *rawReader)
821 // function was called
824 //check if last pad has allready been processed, if not do so
825 if ( fMaxTimeBin>-1 ) ProcessPad();
829 // TVectorF vMeanTime(72);
830 // TVectorF vMeanQ(72);
831 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
832 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
834 //find mean time0 offset for side A and C
835 Double_t time0Side[2]; //time0 for side A:0 and C:1
836 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:1
837 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
838 for ( Int_t iSec = 0; iSec<72; ++iSec ){
839 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
840 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
842 if ( time0SideCount[0] >0 )
843 time0Side[0]/=time0SideCount[0];
844 if ( time0SideCount[1] >0 )
845 time0Side[1]/=time0SideCount[1];
846 // end find time0 offset
848 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
849 for ( Int_t iSec = 0; iSec<72; ++iSec ){
850 //find median and then calculate the mean around it
851 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
852 if ( !hMeanT ) continue;
854 Double_t entries = hMeanT->GetEntries();
856 Short_t *arr = hMeanT->GetArray()+1;
858 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
860 if ( sum>=(entries/2) ) break;
863 Int_t firstBin = fFirstTimeBin+ibin-delta;
864 Int_t lastBin = fFirstTimeBin+ibin+delta;
865 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
866 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
867 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
869 // check boundaries for ebye info of mean time
870 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
871 Int_t vSize=vMeanTime->GetNrows();
872 if ( vSize < fNevents+1 )
873 vMeanTime->ResizeTo(vSize+100);
875 vMeanTime->GetMatrixArray()[fNevents]=median;
878 TVectorF *vTimes = GetPadTimesEvent(iSec);
879 if ( !vTimes ) continue; //continue if no time information for this sector is available
882 AliTPCCalROC calIrocOutliers(0);
883 AliTPCCalROC calOrocOutliers(36);
885 // calculate mean Q of the sector
887 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
888 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
889 if ( vSize < fNevents+1 ) // vSize is the same as for vMeanTime!
890 vMeanQ->ResizeTo(vSize+100);
892 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
894 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
895 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
897 //set values for temporary roc calibration class
899 calIroc->SetValue(iChannel, time);
900 if ( time == 0 ) calIrocOutliers.SetValue(iChannel,1);
903 calOroc->SetValue(iChannel, time);
904 if ( time == 0 ) calOrocOutliers.SetValue(iChannel,1);
907 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
908 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
912 //------------------------------- Debug start ------------------------------
913 if ( fDebugLevel>0 ){
914 if ( !fDebugStreamer ) {
916 TDirectory *backup = gDirectory;
917 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
918 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
925 Float_t q = (*GetPadQEvent(iSec))[iChannel];
926 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
928 UInt_t channel=iChannel;
931 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
932 pad = channel-fROC->GetRowIndexes(sector)[row];
933 padc = pad-(fROC->GetNPads(sector,row)/2);
935 // TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
936 // Form("hSignalD%d.%d.%d",sector,row,pad),
937 // fLastTimeBin-fFirstTimeBin,
938 // fFirstTimeBin,fLastTimeBin);
939 // h1->SetDirectory(0);
941 // for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
942 // h1->Fill(i,fPadSignal(i));
945 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
946 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
947 Double_t t0Side = time0Side[(iSec/18)%2];
948 (*fDebugStreamer) << "DataPad" <<
949 "Event=" << fNevents <<
950 "RunNumber=" << fRunNumber <<
951 "TimeStamp=" << fTimeStamp <<
952 "Sector="<< sector <<
956 "PadSec="<< channel <<
957 "Time0Sec=" << t0Sec <<
958 "Time0Side=" << t0Side <<
969 //----------------------------- Debug end ------------------------------
973 TVectorD paramPol1(3);
974 TVectorD paramPol2(6);
975 TMatrixD matPol1(3,3);
976 TMatrixD matPol2(6,6);
980 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
982 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
983 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
985 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
986 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
989 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
990 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
992 // printf("events: %d -- size: %d\n",fNevents,GetParamArrayPol1(iSec)->GetSize());
994 //------------------------------- Debug start ------------------------------
995 if ( fDebugLevel>0 ){
996 if ( !fDebugStreamer ) {
998 TDirectory *backup = gDirectory;
999 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1000 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1002 (*fDebugStreamer) << "DataRoc" <<
1003 // "Event=" << fEvent <<
1004 "RunNumber=" << fRunNumber <<
1005 "TimeStamp=" << fTimeStamp <<
1007 "hMeanT.=" << hMeanT <<
1008 "median=" << median <<
1009 "paramPol1.=" << ¶mPol1 <<
1010 "paramPol2.=" << ¶mPol2 <<
1011 "matPol1.=" << &matPol1 <<
1012 "matPol2.=" << &matPol2 <<
1013 "chi2Pol1=" << chi2Pol1 <<
1014 "chi2Pol2=" << chi2Pol2 <<
1017 //------------------------------- Debug end ------------------------------
1020 // fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1021 // fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1022 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1023 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1024 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1026 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1027 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1030 fOldRunNumber = fRunNumber;
1035 //_____________________________________________________________________
1036 Bool_t AliTPCCalibCE::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
1039 // Event Processing loop - AliTPCRawStreamFast
1042 Bool_t withInput = kFALSE;
1043 while ( rawStreamFast->NextDDL() ){
1044 while ( rawStreamFast->NextChannel() ){
1045 Int_t isector = rawStreamFast->GetSector(); // current sector
1046 Int_t iRow = rawStreamFast->GetRow(); // current row
1047 Int_t iPad = rawStreamFast->GetPad(); // current pad
1048 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
1049 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
1051 while ( rawStreamFast->NextBunch() ){
1052 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
1053 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
1054 Update(isector,iRow,iPad,iTimeBin+1,signal);
1065 //_____________________________________________________________________
1066 Bool_t AliTPCCalibCE::ProcessEventFast(AliRawReader *rawReader)
1069 // Event processing loop using the fast raw stream algorithm- AliRawReader
1072 //printf("ProcessEventFast - raw reader\n");
1074 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1076 fTimeStamp = eventHeader->Get("Timestamp");
1077 fRunNumber = eventHeader->Get("RunNb");
1079 fEventId = *rawReader->GetEventId();
1081 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader);
1082 Bool_t res=ProcessEventFast(rawStreamFast);
1083 delete rawStreamFast;
1087 //_____________________________________________________________________
1088 Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream)
1091 // Event Processing loop - AliTPCRawStream
1092 // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
1095 rawStream->SetOldRCUFormat(fOldRCUformat);
1099 Bool_t withInput = kFALSE;
1101 while (rawStream->Next()) {
1103 Int_t isector = rawStream->GetSector(); // current sector
1104 Int_t iRow = rawStream->GetRow(); // current row
1105 Int_t iPad = rawStream->GetPad(); // current pad
1106 Int_t iTimeBin = rawStream->GetTime(); // current time bin
1107 Float_t signal = rawStream->GetSignal(); // current ADC signal
1109 Update(isector,iRow,iPad,iTimeBin,signal);
1119 //_____________________________________________________________________
1120 Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
1123 // Event processing loop - AliRawReader
1127 AliTPCRawStream rawStream(rawReader);
1128 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1130 fTimeStamp = eventHeader->Get("Timestamp");
1131 fRunNumber = eventHeader->Get("RunNb");
1133 fEventId = *rawReader->GetEventId();
1136 rawReader->Select("TPC");
1138 return ProcessEvent(&rawStream);
1140 //_____________________________________________________________________
1141 Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
1144 // Event processing loop - date event
1146 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1147 Bool_t result=ProcessEvent(rawReader);
1152 //_____________________________________________________________________
1153 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1154 Int_t nbinsY, Float_t ymin, Float_t ymax,
1155 Char_t *type, Bool_t force)
1158 // return pointer to TH2S histogram of 'type'
1159 // if force is true create a new histogram if it doesn't exist allready
1161 if ( !force || arr->UncheckedAt(sector) )
1162 return (TH2S*)arr->UncheckedAt(sector);
1164 // if we are forced and histogram doesn't exist yet create it
1165 Char_t name[255], title[255];
1167 sprintf(name,"hCalib%s%.2d",type,sector);
1168 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1170 // new histogram with Q calib information. One value for each pad!
1171 TH2S* hist = new TH2S(name,title,
1173 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1174 hist->SetDirectory(0);
1175 arr->AddAt(hist,sector);
1178 //_____________________________________________________________________
1179 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1182 // return pointer to T0 histogram
1183 // if force is true create a new histogram if it doesn't exist allready
1185 TObjArray *arr = &fHistoT0Array;
1186 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1188 //_____________________________________________________________________
1189 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1192 // return pointer to Q histogram
1193 // if force is true create a new histogram if it doesn't exist allready
1195 TObjArray *arr = &fHistoQArray;
1196 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1198 //_____________________________________________________________________
1199 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1202 // return pointer to Q histogram
1203 // if force is true create a new histogram if it doesn't exist allready
1205 TObjArray *arr = &fHistoRMSArray;
1206 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1208 //_____________________________________________________________________
1209 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1210 Char_t *type, Bool_t force)
1213 // return pointer to TH1S histogram
1214 // if force is true create a new histogram if it doesn't exist allready
1216 if ( !force || arr->UncheckedAt(sector) )
1217 return (TH1S*)arr->UncheckedAt(sector);
1219 // if we are forced and histogram doesn't yes exist create it
1220 Char_t name[255], title[255];
1222 sprintf(name,"hCalib%s%.2d",type,sector);
1223 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1225 // new histogram with Q calib information. One value for each pad!
1226 TH1S* hist = new TH1S(name,title,
1227 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1228 hist->SetDirectory(0);
1229 arr->AddAt(hist,sector);
1232 //_____________________________________________________________________
1233 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1236 // return pointer to Q histogram
1237 // if force is true create a new histogram if it doesn't exist allready
1239 TObjArray *arr = &fHistoTmean;
1240 return GetHisto(sector, arr, "LastTmean", force);
1242 //_____________________________________________________________________
1243 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1246 // return pointer to Pad Info from 'arr' for the current event and sector
1247 // if force is true create it if it doesn't exist allready
1249 if ( !force || arr->UncheckedAt(sector) )
1250 return (TVectorF*)arr->UncheckedAt(sector);
1252 TVectorF *vect = new TVectorF(size);
1253 arr->AddAt(vect,sector);
1256 //_____________________________________________________________________
1257 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1260 // return pointer to Pad Times Array for the current event and sector
1261 // if force is true create it if it doesn't exist allready
1263 TObjArray *arr = &fPadTimesArrayEvent;
1264 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1266 //_____________________________________________________________________
1267 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1270 // return pointer to Pad Q Array for the current event and sector
1271 // if force is true create it if it doesn't exist allready
1272 // for debugging purposes only
1275 TObjArray *arr = &fPadQArrayEvent;
1276 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1278 //_____________________________________________________________________
1279 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1282 // return pointer to Pad RMS Array for the current event and sector
1283 // if force is true create it if it doesn't exist allready
1284 // for debugging purposes only
1286 TObjArray *arr = &fPadRMSArrayEvent;
1287 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1289 //_____________________________________________________________________
1290 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1293 // return pointer to Pad RMS Array for the current event and sector
1294 // if force is true create it if it doesn't exist allready
1295 // for debugging purposes only
1297 TObjArray *arr = &fPadPedestalArrayEvent;
1298 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1300 //_____________________________________________________________________
1301 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1304 // return pointer to the EbyE info of the mean arrival time for 'sector'
1305 // if force is true create it if it doesn't exist allready
1307 TObjArray *arr = &fTMeanArrayEvent;
1308 return GetVectSector(sector,arr,100,force);
1310 //_____________________________________________________________________
1311 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1314 // return pointer to the EbyE info of the mean arrival time for 'sector'
1315 // if force is true create it if it doesn't exist allready
1317 TObjArray *arr = &fQMeanArrayEvent;
1318 return GetVectSector(sector,arr,100,force);
1320 //_____________________________________________________________________
1321 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1324 // return pointer to ROC Calibration
1325 // if force is true create a new histogram if it doesn't exist allready
1327 if ( !force || arr->UncheckedAt(sector) )
1328 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1330 // if we are forced and histogram doesn't yes exist create it
1332 // new AliTPCCalROC for T0 information. One value for each pad!
1333 AliTPCCalROC *croc = new AliTPCCalROC(sector);
1334 arr->AddAt(croc,sector);
1337 //_____________________________________________________________________
1338 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1341 // return pointer to Carge ROC Calibration
1342 // if force is true create a new histogram if it doesn't exist allready
1344 TObjArray *arr = &fCalRocArrayT0;
1345 return GetCalRoc(sector, arr, force);
1347 //_____________________________________________________________________
1348 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1351 // return pointer to T0 ROC Calibration
1352 // if force is true create a new histogram if it doesn't exist allready
1354 TObjArray *arr = &fCalRocArrayQ;
1355 return GetCalRoc(sector, arr, force);
1357 //_____________________________________________________________________
1358 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1361 // return pointer to signal width ROC Calibration
1362 // if force is true create a new histogram if it doesn't exist allready
1364 TObjArray *arr = &fCalRocArrayRMS;
1365 return GetCalRoc(sector, arr, force);
1367 //_____________________________________________________________________
1368 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1371 // return pointer to Outliers
1372 // if force is true create a new histogram if it doesn't exist allready
1374 TObjArray *arr = &fCalRocArrayOutliers;
1375 return GetCalRoc(sector, arr, force);
1377 //_____________________________________________________________________
1378 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1381 // return pointer to TObjArray of fit parameters
1382 // if force is true create a new histogram if it doesn't exist allready
1384 if ( !force || arr->UncheckedAt(sector) )
1385 return (TObjArray*)arr->UncheckedAt(sector);
1387 // if we are forced and array doesn't yes exist create it
1389 // new TObjArray for parameters
1390 TObjArray *newArr = new TObjArray;
1391 arr->AddAt(newArr,sector);
1394 //_____________________________________________________________________
1395 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1398 // return pointer to TObjArray of fit parameters from plane fit
1399 // if force is true create a new histogram if it doesn't exist allready
1401 TObjArray *arr = &fParamArrayEventPol1;
1402 return GetParamArray(sector, arr, force);
1404 //_____________________________________________________________________
1405 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1408 // return pointer to TObjArray of fit parameters from parabola fit
1409 // if force is true create a new histogram if it doesn't exist allready
1411 TObjArray *arr = &fParamArrayEventPol2;
1412 return GetParamArray(sector, arr, force);
1414 //_____________________________________________________________________
1415 void AliTPCCalibCE::ResetEvent()
1418 // Reset global counters -- Should be called before each event is processed
1427 fPadTimesArrayEvent.Delete();
1428 fPadQArrayEvent.Delete();
1429 fPadRMSArrayEvent.Delete();
1430 fPadPedestalArrayEvent.Delete();
1432 for ( Int_t i=0; i<72; ++i ){
1433 fVTime0Offset.GetMatrixArray()[i]=0;
1434 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1435 fVMeanQ.GetMatrixArray()[i]=0;
1436 fVMeanQCounter.GetMatrixArray()[i]=0;
1439 //_____________________________________________________________________
1440 void AliTPCCalibCE::ResetPad()
1443 // Reset pad infos -- Should be called after a pad has been processed
1445 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1446 fPadSignal.GetMatrixArray()[i] = 0;
1452 //_____________________________________________________________________
1453 void AliTPCCalibCE::Merge(AliTPCCalibCE *ce)
1456 // Merge ce to the current AliTPCCalibCE
1460 for (Int_t iSec=0; iSec<72; ++iSec){
1461 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1462 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1463 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1467 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1468 TH2S *hRefQ = GetHistoQ(iSec);
1469 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1471 TH2S *hist = new TH2S(*hRefQmerge);
1472 hist->SetDirectory(0);
1473 fHistoQArray.AddAt(hist, iSec);
1475 hRefQmerge->SetDirectory(dir);
1478 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1479 TH2S *hRefT0 = GetHistoT0(iSec);
1480 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1482 TH2S *hist = new TH2S(*hRefT0merge);
1483 hist->SetDirectory(0);
1484 fHistoT0Array.AddAt(hist, iSec);
1486 hRefT0merge->SetDirectory(dir);
1488 if ( hRefRMSmerge ){
1489 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1490 TH2S *hRefRMS = GetHistoRMS(iSec);
1491 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1493 TH2S *hist = new TH2S(*hRefRMSmerge);
1494 hist->SetDirectory(0);
1495 fHistoRMSArray.AddAt(hist, iSec);
1497 hRefRMSmerge->SetDirectory(dir);
1502 // merge time information
1505 Int_t nCEevents = ce->GetNeventsProcessed();
1506 for (Int_t iSec=0; iSec<72; ++iSec){
1507 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1508 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1509 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1510 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1512 TObjArray *arrPol1 = 0x0;
1513 TObjArray *arrPol2 = 0x0;
1514 TVectorF *vMeanTime = 0x0;
1515 TVectorF *vMeanQ = 0x0;
1518 if ( arrPol1CE && arrPol2CE ){
1519 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1520 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1521 arrPol1->Expand(fNevents+nCEevents);
1522 arrPol2->Expand(fNevents+nCEevents);
1524 if ( vMeanTimeCE && vMeanQCE ){
1525 vMeanTime = GetTMeanEvents(iSec);
1526 vMeanQCE = GetQMeanEvents(iSec);
1527 vMeanTime->ResizeTo(fNevents+nCEevents);
1528 vMeanQ->ResizeTo(fNevents+nCEevents);
1532 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1533 if ( arrPol1CE && arrPol2CE ){
1534 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1535 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1536 if ( paramPol1 && paramPol2 ){
1537 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1538 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1541 if ( vMeanTimeCE && vMeanQCE ){
1542 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1543 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1550 TVectorD* eventTimes = ce->GetEventTimes();
1551 TVectorD* eventIds = ce->GetEventIds();
1552 fVEventTime.ResizeTo(fNevents+nCEevents);
1553 fVEventNumber.ResizeTo(fNevents+nCEevents);
1555 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1556 Double_t evTime = eventTimes->GetMatrixArray()[iEvent];
1557 Double_t evId = eventIds->GetMatrixArray()[iEvent];
1558 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1559 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1561 fNevents+=nCEevents; //increase event counter
1564 //_____________________________________________________________________
1565 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1568 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1569 // xVariable: 0-event time, 1-event id, 2-internal event counter
1570 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1571 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1572 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1573 // not used for mean time and mean Q )
1574 // for an example see class description at the beginning
1577 Double_t *x = new Double_t[fNevents];
1578 Double_t *y = new Double_t[fNevents];
1580 TVectorD *xVar = 0x0;
1581 TObjArray *aType = 0x0;
1585 if ( (sector<0) || (sector>71) ) return 0x0;
1586 if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1587 if ( (fitType<0) || (fitType>3) ) return 0x0;
1589 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1590 aType = &fParamArrayEventPol1;
1591 if ( aType->At(sector)==0x0 ) return 0x0;
1593 else if ( fitType==1 ){
1594 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1595 aType = &fParamArrayEventPol2;
1596 if ( aType->At(sector)==0x0 ) return 0x0;
1600 if ( xVariable == 0 ) xVar = &fVEventTime;
1601 if ( xVariable == 1 ) xVar = &fVEventNumber;
1602 if ( xVariable == 2 ) {
1603 xVar = new TVectorD(fNevents);
1604 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1607 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1609 TObjArray *events = (TObjArray*)(aType->At(sector));
1610 if ( events->GetSize()<=ievent ) break;
1611 TVectorD *v = (TVectorD*)(events->At(ievent));
1612 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1613 } else if (fitType == 2) {
1614 Double_t xValue=(*xVar)[ievent];
1615 Double_t yValue=(*GetTMeanEvents(sector))[ievent];
1616 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1617 }else if (fitType == 3) {
1618 Double_t xValue=(*xVar)[ievent];
1619 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1620 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1624 TGraph *gr = new TGraph(npoints);
1625 //sort xVariable increasing
1626 Int_t *sortIndex = new Int_t[npoints];
1627 TMath::Sort(npoints,x,sortIndex);
1628 for (Int_t i=0;i<npoints;++i){
1629 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1633 if ( xVariable == 2 ) delete xVar;
1639 //_____________________________________________________________________
1640 void AliTPCCalibCE::Analyse()
1643 // Calculate calibration constants
1647 TVectorD paramT0(3);
1648 TVectorD paramRMS(3);
1649 TMatrixD dummy(3,3);
1651 for (Int_t iSec=0; iSec<72; ++iSec){
1652 TH2S *hT0 = GetHistoT0(iSec);
1653 if (!hT0 ) continue;
1655 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1656 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1657 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1658 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1660 TH2S *hQ = GetHistoQ(iSec);
1661 TH2S *hRMS = GetHistoRMS(iSec);
1663 Short_t *arrayhQ = hQ->GetArray();
1664 Short_t *arrayhT0 = hT0->GetArray();
1665 Short_t *arrayhRMS = hRMS->GetArray();
1667 UInt_t nChannels = fROC->GetNChannels(iSec);
1675 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1678 Float_t cogTime0 = -1000;
1679 Float_t cogQ = -1000;
1680 Float_t cogRMS = -1000;
1684 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1685 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1686 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1688 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1689 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1690 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1695 //outlier specifications
1696 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1703 rocQ->SetValue(iChannel, cogQ*cogQ);
1704 rocT0->SetValue(iChannel, cogTime0);
1705 rocRMS->SetValue(iChannel, cogRMS);
1706 rocOut->SetValue(iChannel, cogOut);
1710 if ( fDebugLevel > 0 ){
1711 if ( !fDebugStreamer ) {
1713 TDirectory *backup = gDirectory;
1714 fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1715 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1718 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1719 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1720 padc = pad-(fROC->GetNPads(iSec,row)/2);
1722 (*fDebugStreamer) << "DataEnd" <<
1723 "Sector=" << iSec <<
1727 "PadSec=" << iChannel <<
1729 "T0=" << cogTime0 <<
1738 if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
1739 // delete fDebugStreamer;
1740 // fDebugStreamer = 0x0;
1742 //_____________________________________________________________________
1743 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1746 // Write class to file
1755 option = "recreate";
1757 TDirectory *backup = gDirectory;
1758 TFile f(filename,option.Data());
1760 if ( !sDir.IsNull() ){
1761 f.mkdir(sDir.Data());
1767 if ( backup ) backup->cd();