]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCalibCE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
calibration updates (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibCE.cxx
CommitLineData
75d8233f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
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 **************************************************************************/
15
75d8233f 16/* $Id$ */
17
3cd27a08 18////////////////////////////////////////////////////////////////////////////////////////
19// //
20// Implementation of the TPC Central Electrode calibration //
21// //
22// Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
23// //
24////////////////////////////////////////////////////////////////////////////////////////
7fb602b1 25//
26//
27// *************************************************************************************
28// * Class Description *
29// *************************************************************************************
30//
31/* BEGIN_HTML
32 <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
33 using laser runs.</h4>
34
35 The information retrieved is
36 <ul style="list-style-type: square;">
37 <li>Time arrival from the CE</li>
38 <li>Signal width</li>
39 <li>Signal sum</li>
40 </ul>
41
42<h4>Overview:</h4>
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>
47 </ol>
48
49 <h3><a name="working">I. Working principle</a></h3>
50
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>
53
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
59 stored in fPadSignal.
60 </li>
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>
71 </ul>
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>
76 </ul>
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>
83 </ul>
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>
86 </ol>
87 </ul>
88 </ul>
89
90 <h4>At the end of each event the EndEvent() function is called</h4>
91
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>
100 </ul>
101 </ul>
102
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>
111 </ul>
112 </ul>
113
114 <h3><a name="user">II. User interface for filling data</a></h3>
115
116 <h4>To Fill information one of the following functions can be used:</h4>
117
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>
123 </ul>
124 <br />
125
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>
130 </ul>
131 <br />
132
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>
137 </ul>
138 <br />
139
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>
145 </ul>
146 </ul>
147
148 <h4>It is also possible to merge two independently taken calibrations using the function</h4>
149
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>
156 </ul>
157 </ul>
158
159
160 <h4>example: filling data using root raw data:</h4>
161 <pre>
162 void fillCE(Char_t *filename)
163 {
164 rawReader = new AliRawReaderRoot(fileName);
165 if ( !rawReader ) return;
166 AliTPCCalibCE *calib = new AliTPCCalibCE;
167 while (rawReader->NextEvent()){
168 calib->ProcessEvent(rawReader);
169 }
170 calib->Analyse();
171 calib->DumpToFile("CEData.root");
172 delete rawReader;
173 delete calib;
174 }
175 </pre>
176
177 <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
178
179 <h4><a name="info:stored">III.1 Stored information</a></h4>
180 <ul style="list-style-type: none;">
181 <li>Histograms:</li>
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>
186 </ul>
187 <br />
188
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>
195 </ul>
196 <br />
197
198 <li>For each event the following information is stored:</li>
199
200 <ul style="list-style-type: square;">
201 <li>event time ( TVectorD fVEventTime )</li>
202 <li>event id ( TVectorD fVEventNumber )</li>
203 <br />
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>
208 </ul>
209 </ul>
210
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>
218 </ul>
219 <br />
220
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>
226 </ul>
227 <br />
228
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>
234 </ul>
235 </ul>
236
237 <h4>example for visualisation:</h4>
238 <pre>
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");
244
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
250
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");
256 </pre>
257END_HTML */
258//////////////////////////////////////////////////////////////////////////////////////
259
260
3cd27a08 261//Root includes
262#include <TObjArray.h>
263#include <TH1F.h>
264#include <TH2S.h>
265#include <TString.h>
266#include <TVectorF.h>
267#include <TVectorD.h>
268#include <TMatrixD.h>
269#include <TMath.h>
270#include <TGraph.h>
271#include <TString.h>
ac940b58 272#include <TMap.h>
3cd27a08 273#include <TDirectory.h>
274#include <TSystem.h>
275#include <TFile.h>
276
277//AliRoot includes
4c6d06dc 278#include "AliLog.h"
3cd27a08 279#include "AliRawReader.h"
280#include "AliRawReaderRoot.h"
281#include "AliRawReaderDate.h"
282#include "AliRawEventHeaderBase.h"
283#include "AliTPCRawStream.h"
08205ed7 284#include "AliTPCRawStreamFast.h"
3cd27a08 285#include "AliTPCcalibDB.h"
286#include "AliTPCCalROC.h"
287#include "AliTPCCalPad.h"
288#include "AliTPCROC.h"
289#include "AliTPCParam.h"
290#include "AliTPCCalibCE.h"
291#include "AliMathBase.h"
292#include "TTreeStream.h"
293
294//date
295#include "event.h"
296ClassImp(AliTPCCalibCE)
297
298
7fb602b1 299AliTPCCalibCE::AliTPCCalibCE() :
75d8233f 300 TObject(),
301 fFirstTimeBin(650),
302 fLastTimeBin(1000),
303 fNbinsT0(200),
304 fXminT0(-5),
305 fXmaxT0(5),
306 fNbinsQ(200),
bf57d87d 307 fXminQ(1),
308 fXmaxQ(40),
75d8233f 309 fNbinsRMS(100),
bf57d87d 310 fXminRMS(0.1),
311 fXmaxRMS(5.1),
4c6d06dc 312 fPeakMinus(2),
313 fPeakPlus(3),
314 fNoiseThresholdMax(5.),
315 fNoiseThresholdSum(8.),
316 fIsZeroSuppressed(kFALSE),
75d8233f 317 fLastSector(-1),
930a9868 318 fSecRejectRatio(.4),
75d8233f 319 fROC(AliTPCROC::Instance()),
04049c81 320 fMapping(NULL),
75d8233f 321 fParam(new AliTPCParam),
322 fPedestalTPC(0x0),
323 fPadNoiseTPC(0x0),
324 fPedestalROC(0x0),
325 fPadNoiseROC(0x0),
75d8233f 326 fCalRocArrayT0(72),
ef7f7670 327 fCalRocArrayT0Err(72),
75d8233f 328 fCalRocArrayQ(72),
329 fCalRocArrayRMS(72),
330 fCalRocArrayOutliers(72),
331 fHistoQArray(72),
332 fHistoT0Array(72),
333 fHistoRMSArray(72),
ef7f7670 334 fMeanT0rms(0),
335 fMeanQrms(0),
336 fMeanRMSrms(0),
75d8233f 337 fHistoTmean(72),
75d8233f 338 fParamArrayEventPol1(72),
339 fParamArrayEventPol2(72),
7fb602b1 340 fTMeanArrayEvent(72),
341 fQMeanArrayEvent(72),
342 fVEventTime(10),
343 fVEventNumber(10),
75d8233f 344 fNevents(0),
345 fTimeStamp(0),
7fb602b1 346 fEventId(-1),
75d8233f 347 fRunNumber(-1),
348 fOldRunNumber(-1),
349 fPadTimesArrayEvent(72),
350 fPadQArrayEvent(72),
351 fPadRMSArrayEvent(72),
352 fPadPedestalArrayEvent(72),
353 fCurrentChannel(-1),
354 fCurrentSector(-1),
355 fCurrentRow(-1),
356 fMaxPadSignal(-1),
357 fMaxTimeBin(-1),
358 fPadSignal(1024),
359 fPadPedestal(0),
360 fPadNoise(0),
361 fVTime0Offset(72),
362 fVTime0OffsetCounter(72),
bf57d87d 363 fVMeanQ(72),
364 fVMeanQCounter(72),
ac940b58 365 fCurrentCETimeRef(0),
3cd27a08 366// fEvent(-1),
75d8233f 367 fDebugStreamer(0x0),
368 fDebugLevel(0)
369{
370 //
371 // AliTPCSignal default constructor
372 //
373// fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
4c6d06dc 374 fParam->Update();
75d8233f 375}
376//_____________________________________________________________________
377AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
378 TObject(sig),
379 fFirstTimeBin(sig.fFirstTimeBin),
380 fLastTimeBin(sig.fLastTimeBin),
381 fNbinsT0(sig.fNbinsT0),
382 fXminT0(sig.fXminT0),
383 fXmaxT0(sig.fXmaxT0),
384 fNbinsQ(sig.fNbinsQ),
385 fXminQ(sig.fXminQ),
386 fXmaxQ(sig.fXmaxQ),
387 fNbinsRMS(sig.fNbinsRMS),
388 fXminRMS(sig.fXminRMS),
389 fXmaxRMS(sig.fXmaxRMS),
4c6d06dc 390 fPeakMinus(sig.fPeakMinus),
391 fPeakPlus(sig.fPeakPlus),
392 fNoiseThresholdMax(sig.fNoiseThresholdMax),
393 fNoiseThresholdSum(sig.fNoiseThresholdSum),
394 fIsZeroSuppressed(sig.fIsZeroSuppressed),
75d8233f 395 fLastSector(-1),
930a9868 396 fSecRejectRatio(.4),
75d8233f 397 fROC(AliTPCROC::Instance()),
04049c81 398 fMapping(NULL),
75d8233f 399 fParam(new AliTPCParam),
400 fPedestalTPC(0x0),
401 fPadNoiseTPC(0x0),
402 fPedestalROC(0x0),
403 fPadNoiseROC(0x0),
75d8233f 404 fCalRocArrayT0(72),
ef7f7670 405 fCalRocArrayT0Err(72),
75d8233f 406 fCalRocArrayQ(72),
407 fCalRocArrayRMS(72),
408 fCalRocArrayOutliers(72),
409 fHistoQArray(72),
410 fHistoT0Array(72),
411 fHistoRMSArray(72),
ef7f7670 412 fMeanT0rms(sig.fMeanT0rms),
413 fMeanQrms(sig.fMeanQrms),
414 fMeanRMSrms(sig.fMeanRMSrms),
75d8233f 415 fHistoTmean(72),
75d8233f 416 fParamArrayEventPol1(72),
417 fParamArrayEventPol2(72),
7fb602b1 418 fTMeanArrayEvent(72),
419 fQMeanArrayEvent(72),
75d8233f 420 fVEventTime(1000),
421 fVEventNumber(1000),
422 fNevents(sig.fNevents),
423 fTimeStamp(0),
7fb602b1 424 fEventId(-1),
75d8233f 425 fRunNumber(-1),
426 fOldRunNumber(-1),
427 fPadTimesArrayEvent(72),
428 fPadQArrayEvent(72),
429 fPadRMSArrayEvent(72),
430 fPadPedestalArrayEvent(72),
431 fCurrentChannel(-1),
432 fCurrentSector(-1),
433 fCurrentRow(-1),
434 fMaxPadSignal(-1),
435 fMaxTimeBin(-1),
436 fPadSignal(1024),
437 fPadPedestal(0),
bf57d87d 438 fPadNoise(0),
75d8233f 439 fVTime0Offset(72),
440 fVTime0OffsetCounter(72),
bf57d87d 441 fVMeanQ(72),
442 fVMeanQCounter(72),
ac940b58 443 fCurrentCETimeRef(0),
3cd27a08 444// fEvent(-1),
75d8233f 445 fDebugStreamer(0x0),
446 fDebugLevel(sig.fDebugLevel)
447{
ac940b58 448 //
7fb602b1 449 // AliTPCSignal copy constructor
ac940b58 450 //
75d8233f 451
ac940b58 452 for (Int_t iSec = 0; iSec < 72; ++iSec){
453 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
454 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
455 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
456 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
75d8233f 457
ac940b58 458 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
459 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
460 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
75d8233f 461
ac940b58 462 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
463 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
464 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
465 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
75d8233f 466
ac940b58 467 if ( hQ != 0x0 ){
7fb602b1 468// TDirectory *dir = hQ->GetDirectory();
469// hQ->SetDirectory(0);
ac940b58 470 TH2S *hNew = new TH2S(*hQ);
471 hNew->SetDirectory(0);
472 fHistoQArray.AddAt(hNew,iSec);
7fb602b1 473// hQ->SetDirectory(dir);
ac940b58 474 }
475 if ( hT0 != 0x0 ){
7fb602b1 476// TDirectory *dir = hT0->GetDirectory();
477// hT0->SetDirectory(0);
ac940b58 478 TH2S *hNew = new TH2S(*hT0);
479 hNew->SetDirectory(0);
480 fHistoT0Array.AddAt(hNew,iSec);
7fb602b1 481// hT0->SetDirectory(dir);
ac940b58 482 }
483 if ( hRMS != 0x0 ){
7fb602b1 484// TDirectory *dir = hRMS->GetDirectory();
485// hRMS->SetDirectory(0);
ac940b58 486 TH2S *hNew = new TH2S(*hRMS);
487 hNew->SetDirectory(0);
488 fHistoRMSArray.AddAt(hNew,iSec);
7fb602b1 489// hRMS->SetDirectory(dir);
75d8233f 490 }
ac940b58 491 }
75d8233f 492
7fb602b1 493 //copy fit parameters event by event
ac940b58 494 TObjArray *arr=0x0;
495 for (Int_t iSec=0; iSec<72; ++iSec){
496 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
497 if ( arr ){
498 TObjArray *arrEvents = new TObjArray(arr->GetSize());
499 fParamArrayEventPol1.AddAt(arrEvents, iSec);
500 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
501 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
502 arrEvents->AddAt(new TVectorD(*vec),iEvent);
503 }
7fb602b1 504
ac940b58 505 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
506 if ( arr ){
507 TObjArray *arrEvents = new TObjArray(arr->GetSize());
508 fParamArrayEventPol2.AddAt(arrEvents, iSec);
509 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
510 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
511 arrEvents->AddAt(new TVectorD(*vec),iEvent);
7fb602b1 512 }
513
ac940b58 514 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
515 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
516 if ( vMeanTime )
517 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
518 if ( vMeanQ )
519 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
520 }
521
7fb602b1 522
ac940b58 523 fVEventTime.ResizeTo(sig.fVEventTime);
524 fVEventNumber.ResizeTo(sig.fVEventNumber);
525 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
526 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
7fb602b1 527
ac940b58 528 fParam->Update();
75d8233f 529}
ac940b58 530//_____________________________________________________________________
531AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
532 TObject(),
533 fFirstTimeBin(650),
534 fLastTimeBin(1000),
535 fNbinsT0(200),
536 fXminT0(-5),
537 fXmaxT0(5),
538 fNbinsQ(200),
539 fXminQ(1),
540 fXmaxQ(40),
541 fNbinsRMS(100),
542 fXminRMS(0.1),
543 fXmaxRMS(5.1),
544 fPeakMinus(2),
545 fPeakPlus(3),
546 fNoiseThresholdMax(5.),
547 fNoiseThresholdSum(8.),
548 fIsZeroSuppressed(kFALSE),
549 fLastSector(-1),
550 fSecRejectRatio(.4),
551 fROC(AliTPCROC::Instance()),
552 fMapping(NULL),
553 fParam(new AliTPCParam),
554 fPedestalTPC(0x0),
555 fPadNoiseTPC(0x0),
556 fPedestalROC(0x0),
557 fPadNoiseROC(0x0),
558 fCalRocArrayT0(72),
559 fCalRocArrayT0Err(72),
560 fCalRocArrayQ(72),
561 fCalRocArrayRMS(72),
562 fCalRocArrayOutliers(72),
563 fHistoQArray(72),
564 fHistoT0Array(72),
565 fHistoRMSArray(72),
566 fMeanT0rms(0),
567 fMeanQrms(0),
568 fMeanRMSrms(0),
569 fHistoTmean(72),
570 fParamArrayEventPol1(72),
571 fParamArrayEventPol2(72),
572 fTMeanArrayEvent(72),
573 fQMeanArrayEvent(72),
574 fVEventTime(10),
575 fVEventNumber(10),
576 fNevents(0),
577 fTimeStamp(0),
578 fEventId(-1),
579 fRunNumber(-1),
580 fOldRunNumber(-1),
581 fPadTimesArrayEvent(72),
582 fPadQArrayEvent(72),
583 fPadRMSArrayEvent(72),
584 fPadPedestalArrayEvent(72),
585 fCurrentChannel(-1),
586 fCurrentSector(-1),
587 fCurrentRow(-1),
588 fMaxPadSignal(-1),
589 fMaxTimeBin(-1),
590 fPadSignal(1024),
591 fPadPedestal(0),
592 fPadNoise(0),
593 fVTime0Offset(72),
594 fVTime0OffsetCounter(72),
595 fVMeanQ(72),
596 fVMeanQCounter(72),
597 fCurrentCETimeRef(0),
598 // fEvent(-1),
599 fDebugStreamer(0x0),
600 fDebugLevel(0)
601{
602 //
603 // constructor which uses a tmap as input to set some specific parameters
604 //
605 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
606 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
607 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
608 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
609 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
610 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
611 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
612 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
613 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
614 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
615 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
616 if (config->GetValue("PeakMinus")) fPeakMinus = ((TObjString*)config->GetValue("PeakMinus"))->GetString().Atoi();
617 if (config->GetValue("PeakPlus")) fPeakPlus = ((TObjString*)config->GetValue("PeakPlus"))->GetString().Atoi();
618 if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
619 if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
620 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
621 if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
622
623 fParam->Update();
624}
625
75d8233f 626//_____________________________________________________________________
627AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
628{
629 //
630 // assignment operator
631 //
632 if (&source == this) return *this;
633 new (this) AliTPCCalibCE(source);
634
635 return *this;
636}
637//_____________________________________________________________________
638AliTPCCalibCE::~AliTPCCalibCE()
639{
640 //
641 // destructor
642 //
643
644 fCalRocArrayT0.Delete();
ef7f7670 645 fCalRocArrayT0Err.Delete();
75d8233f 646 fCalRocArrayQ.Delete();
647 fCalRocArrayRMS.Delete();
7fb602b1 648 fCalRocArrayOutliers.Delete();
75d8233f 649
650 fHistoQArray.Delete();
651 fHistoT0Array.Delete();
652 fHistoRMSArray.Delete();
653
7fb602b1 654 fHistoTmean.Delete();
655
656 fParamArrayEventPol1.Delete();
657 fParamArrayEventPol2.Delete();
658 fTMeanArrayEvent.Delete();
659 fQMeanArrayEvent.Delete();
660
75d8233f 661 fPadTimesArrayEvent.Delete();
662 fPadQArrayEvent.Delete();
663 fPadRMSArrayEvent.Delete();
664 fPadPedestalArrayEvent.Delete();
665
666 if ( fDebugStreamer) delete fDebugStreamer;
667// if ( fHTime0 ) delete fHTime0;
d8faeea7 668// delete fROC;
75d8233f 669 delete fParam;
670}
671//_____________________________________________________________________
7fb602b1 672Int_t AliTPCCalibCE::Update(const Int_t icsector,
75d8233f 673 const Int_t icRow,
674 const Int_t icPad,
675 const Int_t icTimeBin,
676 const Float_t csignal)
677{
ac940b58 678 //
679 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
680 // no extra analysis necessary. Assumes knowledge of the signal shape!
681 // assumes that it is looped over consecutive time bins of one pad
682 //
4c6d06dc 683
ac940b58 684 //temp
4c6d06dc 685
b401648b 686 if (icRow<0) return 0;
687 if (icPad<0) return 0;
688 if (icTimeBin<0) return 0;
ac940b58 689 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
75d8233f 690
ac940b58 691 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
75d8233f 692
ac940b58 693 //init first pad and sector in this event
694 if ( fCurrentChannel == -1 ) {
695 fLastSector=-1;
696 fCurrentChannel = iChannel;
697 fCurrentSector = icsector;
698 fCurrentRow = icRow;
699 }
75d8233f 700
ac940b58 701 //process last pad if we change to a new one
702 if ( iChannel != fCurrentChannel ){
703 ProcessPad();
704 fLastSector=fCurrentSector;
705 fCurrentChannel = iChannel;
706 fCurrentSector = icsector;
707 fCurrentRow = icRow;
708 }
75d8233f 709
ac940b58 710 //fill signals for current pad
711 fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
712 if ( csignal > fMaxPadSignal ){
713 fMaxPadSignal = csignal;
714 fMaxTimeBin = icTimeBin;
715 }
716 return 0;
75d8233f 717}
718//_____________________________________________________________________
719void AliTPCCalibCE::FindPedestal(Float_t part)
720{
ac940b58 721 //
75d8233f 722 // find pedestal and noise for the current pad. Use either database or
723 // truncated mean with part*100%
ac940b58 724 //
725 Bool_t noPedestal = kTRUE;
7fb602b1 726
727 //use pedestal database if set
ac940b58 728 if (fPedestalTPC&&fPadNoiseTPC){
75d8233f 729 //only load new pedestals if the sector has changed
ac940b58 730 if ( fCurrentSector!=fLastSector ){
731 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
732 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
733 }
75d8233f 734
ac940b58 735 if ( fPedestalROC&&fPadNoiseROC ){
736 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
737 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
738 noPedestal = kFALSE;
75d8233f 739 }
740
ac940b58 741 }
742
75d8233f 743 //if we are not running with pedestal database, or for the current sector there is no information
744 //available, calculate the pedestal and noise on the fly
ac940b58 745 if ( noPedestal ) {
746 fPadPedestal = 0;
747 fPadNoise = 0;
748 if ( fIsZeroSuppressed ) return;
749 const Int_t kPedMax = 100; //maximum pedestal value
750 Float_t max = 0;
751 Float_t maxPos = 0;
752 Int_t median = -1;
753 Int_t count0 = 0;
754 Int_t count1 = 0;
755 //
756 Float_t padSignal=0;
757 //
758 UShort_t histo[kPedMax];
759 memset(histo,0,kPedMax*sizeof(UShort_t));
bf57d87d 760
7fb602b1 761 //fill pedestal histogram
ac940b58 762 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
763 padSignal = fPadSignal.GetMatrixArray()[i];
764 if (padSignal<=0) continue;
765 if (padSignal>max && i>10) {
766 max = padSignal;
767 maxPos = i;
768 }
769 if (padSignal>kPedMax-1) continue;
770 histo[int(padSignal+0.5)]++;
771 count0++;
772 }
7fb602b1 773 //find median
ac940b58 774 for (Int_t i=1; i<kPedMax; ++i){
775 if (count1<count0*0.5) median=i;
776 count1+=histo[i];
777 }
75d8233f 778 // truncated mean
ac940b58 779 //
780 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
781 //
782 for (Int_t idelta=1; idelta<10; ++idelta){
783 if (median-idelta<=0) continue;
784 if (median+idelta>kPedMax) continue;
785 if (count<part*count1){
786 count+=histo[median-idelta];
787 mean +=histo[median-idelta]*(median-idelta);
788 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
789 count+=histo[median+idelta];
790 mean +=histo[median+idelta]*(median+idelta);
791 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
792 }
793 }
794 if ( count > 0 ) {
795 mean/=count;
796 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
797 fPadPedestal = mean;
798 fPadNoise = rms;
75d8233f 799 }
ac940b58 800 }
75d8233f 801}
802//_____________________________________________________________________
ac940b58 803void AliTPCCalibCE::UpdateCETimeRef()
804{
805 // Find the time reference of the last valid CE signal in sector
806 // for irocs of the A-Side the reference of the corresponging OROC is returned
807 // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
808 if ( fLastSector == fCurrentSector ) return;
809 Int_t sector=fCurrentSector;
810 if ( sector < 18 ) sector+=36;
811 fCurrentCETimeRef=0;
812 TVectorF *vtRef = GetTMeanEvents(sector);
813 if ( !vtRef ) return;
814 Int_t vtRefSize= vtRef->GetNrows();
815 if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
816 else vtRefSize=fNevents;
817 while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
818 fCurrentCETimeRef=(*vtRef)[vtRefSize];
819 AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
820}
821//_____________________________________________________________________
75d8233f 822void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
823{
ac940b58 824 //
75d8233f 825 // Find position, signal width and height of the CE signal (last signal)
826 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
827 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
ac940b58 828 //
75d8233f 829
ac940b58 830 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
831 Int_t cemaxpos = 0;
832 Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise; // threshold for the signal sum
833 const Int_t kCemin = 4; // range for the analysis of the ce signal +- channels from the peak
834 const Int_t kCemax = 7;
75d8233f 835
ac940b58 836 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
bf57d87d 837
75d8233f 838 // find maximum closest to the sector mean from the last event
ac940b58 839 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
7fb602b1 840 // get sector mean of last event
ac940b58 841 Float_t tmean = fCurrentCETimeRef;
842 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
843 minDist = tmean-maxima[imax];
844 cemaxpos = (Int_t)maxima[imax];
75d8233f 845 }
ac940b58 846 }
75d8233f 847
ac940b58 848 if (cemaxpos!=0){
849 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
850 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
851 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
852 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
853 if (signal>0) {
854 ceTime+=signal*(i+0.5);
855 ceRMS +=signal*(i+0.5)*(i+0.5);
856 ceQsum+=signal;
857 }
858 }
75d8233f 859 }
ac940b58 860 }
861 if (ceQmax&&ceQsum>ceSumThreshold) {
862 ceTime/=ceQsum;
863 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
864 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
865 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
866
867 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
868 // the pick-up signal should scale with the pad area. In addition
869 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
870 // ratio 2/3. The pad area we express in cm2. We normalise the signal
871 // to the OROC signal (factor 2/3 for the IROCs).
872 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
873 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
874
875 ceQsum/=norm;
876 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
877 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
878 } else {
879 ceQmax=0;
880 ceTime=0;
881 ceRMS =0;
882 ceQsum=0;
883 }
884 param[0] = ceQmax;
885 param[1] = ceTime;
886 param[2] = ceRMS;
887 qSum = ceQsum;
75d8233f 888}
889//_____________________________________________________________________
3cd27a08 890Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
75d8233f 891{
ac940b58 892 //
75d8233f 893 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
894 // and 'tplus' timebins after 'pos'
ac940b58 895 //
896 if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
897 for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
898 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
899 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
900 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
901 ++iTime2;
902 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
903 }
904 return kTRUE;
75d8233f 905}
906//_____________________________________________________________________
907void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
908{
ac940b58 909 //
75d8233f 910 // Find local maxima on the pad signal and Histogram them
ac940b58 911 //
4c6d06dc 912 Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal
ac940b58 913 Int_t count = 0;
4c6d06dc 914// Int_t tminus = 2;
915// Int_t tplus = 3;
ac940b58 916 for (Int_t i=fLastTimeBin-fPeakPlus-1; i>=fFirstTimeBin+fPeakMinus; --i){
917 if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,fPeakMinus,fPeakPlus) ){
918 if (count<maxima.GetNrows()){
919 maxima.GetMatrixArray()[count++]=i;
920 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
921 }
75d8233f 922 }
ac940b58 923 }
75d8233f 924}
925//_____________________________________________________________________
7fb602b1 926void AliTPCCalibCE::ProcessPad()
75d8233f 927{
928 //
929 // Process data of current pad
930 //
931 FindPedestal();
932
7fb602b1 933 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
934 // + central electrode and possibly post peaks from the CE signal
935 // however if we are on a high noise pad a lot more peaks due to the noise might occur
75d8233f 936 FindLocalMaxima(maxima);
937 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
938
ac940b58 939 UpdateCETimeRef(); // update the time refenrence for the current sector
940 if ( fCurrentCETimeRef==0 ) return; //return if we don't have time 0 info, eg if only one side has laser
75d8233f 941 TVectorD param(3);
3cd27a08 942 Float_t qSum;
943 FindCESignal(param, qSum, maxima);
75d8233f 944
945 Double_t meanT = param[1];
946 Double_t sigmaT = param[2];
947
948 //Fill Event T0 counter
bf57d87d 949 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
75d8233f 950
951 //Fill Q histogram
3cd27a08 952 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
75d8233f 953
954 //Fill RMS histogram
955 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
956
957
958 //Fill debugging info
959 if ( fDebugLevel>0 ){
bf57d87d 960 (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
961 (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
3cd27a08 962 (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
75d8233f 963 }
964
965 ResetPad();
966}
967//_____________________________________________________________________
7fb602b1 968void AliTPCCalibCE::EndEvent()
75d8233f 969{
ac940b58 970 // Process data of current pad
971 // The Functions 'SetTimeStamp' and 'SetRunNumber' should be called
972 // before the EndEvent function to set the event timestamp and number!!!
973 // This is automatically done if the ProcessEvent(AliRawReader *rawReader)
974 // function was called
75d8233f 975
ac940b58 976 //check if last pad has allready been processed, if not do so
977 if ( fMaxTimeBin>-1 ) ProcessPad();
75d8233f 978
ac940b58 979 AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
4c6d06dc 980
ac940b58 981 TVectorD param(3);
982 TMatrixD dummy(3,3);
7fb602b1 983// TVectorF vMeanTime(72);
984// TVectorF vMeanQ(72);
ac940b58 985 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
986 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
987
988 //find mean time0 offset for side A and C
989 //use only orocs due to the better statistics
990 Double_t time0Side[2]; //time0 for side A:0 and C:1
991 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:1
992 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
993 for ( Int_t iSec = 36; iSec<72; ++iSec ){
994 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
995 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
996 }
997 if ( time0SideCount[0] >0 )
998 time0Side[0]/=time0SideCount[0];
999 if ( time0SideCount[1] >0 )
1000 time0Side[1]/=time0SideCount[1];
75d8233f 1001 // end find time0 offset
ac940b58 1002 AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1003 Int_t nSecMeanT=0;
1004 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1005 for ( Int_t iSec = 0; iSec<72; ++iSec ){
1006 AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1007 //find median and then calculate the mean around it
1008 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
1009 if ( !hMeanT ) continue;
1010 //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1011 if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1012 hMeanT->Reset();
1013 AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1014 continue;
1015 }
1016
1017 Double_t entries = hMeanT->GetEffectiveEntries();
1018 Double_t sum = 0;
1019 Short_t *arr = hMeanT->GetArray()+1;
1020 Int_t ibin=0;
1021 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1022 sum+=arr[ibin];
1023 if ( sum>=(entries/2.) ) break;
1024 }
1025 Int_t delta = 4;
1026 Int_t firstBin = fFirstTimeBin+ibin-delta;
1027 Int_t lastBin = fFirstTimeBin+ibin+delta;
1028 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1029 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
1030 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1031
7fb602b1 1032 // check boundaries for ebye info of mean time
ac940b58 1033 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1034 Int_t vSize=vMeanTime->GetNrows();
1035 if ( vSize < fNevents+1 ){
1036 vMeanTime->ResizeTo(vSize+100);
1037 }
1038
1039 vMeanTime->GetMatrixArray()[fNevents]=median;
1040 nSecMeanT++;
1041 // end find median
1042
1043 TVectorF *vTimes = GetPadTimesEvent(iSec);
1044 if ( !vTimes ) continue; //continue if no time information for this sector is available
1045
1046 AliTPCCalROC calIrocOutliers(0);
1047 AliTPCCalROC calOrocOutliers(36);
1048
1049 // calculate mean Q of the sector
1050 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1051 vSize=vMeanQ->GetNrows();
1052 if ( vSize < fNevents+1 ){
1053 vMeanQ->ResizeTo(vSize+100);
1054 }
1055 Float_t meanQ = 0;
1056 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1057 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1058
1059 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1060 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
75d8233f 1061
1062 //set values for temporary roc calibration class
ac940b58 1063 if ( iSec < 36 ) {
1064 calIroc->SetValue(iChannel, time);
1065 if ( time == 0 ) calIrocOutliers.SetValue(iChannel,1);
bf57d87d 1066
ac940b58 1067 } else {
1068 calOroc->SetValue(iChannel, time);
1069 if ( time == 0 ) calOrocOutliers.SetValue(iChannel,1);
1070 }
75d8233f 1071
ac940b58 1072 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1073 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
75d8233f 1074
1075
1076
1077 //------------------------------- Debug start ------------------------------
ac940b58 1078 if ( fDebugLevel>0 ){
1079 if ( !fDebugStreamer ) {
75d8233f 1080 //debug stream
ac940b58 1081 TDirectory *backup = gDirectory;
1082 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1083 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1084 }
75d8233f 1085
ac940b58 1086 Int_t row=0;
1087 Int_t pad=0;
1088 Int_t padc=0;
75d8233f 1089
ac940b58 1090 Float_t q = (*GetPadQEvent(iSec))[iChannel];
1091 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
75d8233f 1092
ac940b58 1093 UInt_t channel=iChannel;
1094 Int_t sector=iSec;
75d8233f 1095
ac940b58 1096 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1097 pad = channel-fROC->GetRowIndexes(sector)[row];
1098 padc = pad-(fROC->GetNPads(sector,row)/2);
75d8233f 1099
1100// TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1101// Form("hSignalD%d.%d.%d",sector,row,pad),
1102// fLastTimeBin-fFirstTimeBin,
1103// fFirstTimeBin,fLastTimeBin);
1104// h1->SetDirectory(0);
ac940b58 1105 //
7fb602b1 1106// for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
75d8233f 1107// h1->Fill(i,fPadSignal(i));
1108
ac940b58 1109 Double_t t0Sec = 0;
1110 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1111 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1112 Double_t t0Side = time0Side[(iSec/18)%2];
1113 (*fDebugStreamer) << "DataPad" <<
1114 "Event=" << fNevents <<
1115 "RunNumber=" << fRunNumber <<
1116 "TimeStamp=" << fTimeStamp <<
1117 "Sector="<< sector <<
1118 "Row=" << row<<
1119 "Pad=" << pad <<
1120 "PadC=" << padc <<
1121 "PadSec="<< channel <<
1122 "Time0Sec=" << t0Sec <<
1123 "Time0Side=" << t0Side <<
1124 "Time=" << time <<
1125 "RMS=" << rms <<
1126 "Sum=" << q <<
1127 "MeanQ=" << meanQ <<
bf57d87d 1128 // "hist.=" << h1 <<
ac940b58 1129 "\n";
75d8233f 1130
bf57d87d 1131 // delete h1;
1132
ac940b58 1133 }
75d8233f 1134 //----------------------------- Debug end ------------------------------
ac940b58 1135 }// end channel loop
1136
1137 TVectorD paramPol1(3);
1138 TVectorD paramPol2(6);
1139 TMatrixD matPol1(3,3);
1140 TMatrixD matPol2(6,6);
1141 Float_t chi2Pol1=0;
1142 Float_t chi2Pol2=0;
1143
1144 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1145 if ( iSec < 36 ){
1146 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1147 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1148 } else {
1149 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1150 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1151 }
75d8233f 1152
ac940b58 1153 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1154 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1155 }
75d8233f 1156
1157 //------------------------------- Debug start ------------------------------
ac940b58 1158 if ( fDebugLevel>0 ){
1159 if ( !fDebugStreamer ) {
bf57d87d 1160 //debug stream
ac940b58 1161 TDirectory *backup = gDirectory;
1162 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1163 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1164 }
1165 (*fDebugStreamer) << "DataRoc" <<
3cd27a08 1166// "Event=" << fEvent <<
ac940b58 1167 "RunNumber=" << fRunNumber <<
1168 "TimeStamp=" << fTimeStamp <<
1169 "Sector="<< iSec <<
1170 "hMeanT.=" << hMeanT <<
1171 "median=" << median <<
1172 "paramPol1.=" << &paramPol1 <<
1173 "paramPol2.=" << &paramPol2 <<
1174 "matPol1.=" << &matPol1 <<
1175 "matPol2.=" << &matPol2 <<
1176 "chi2Pol1=" << chi2Pol1 <<
1177 "chi2Pol2=" << chi2Pol2 <<
1178 "\n";
1179 }
75d8233f 1180 //------------------------------- Debug end ------------------------------
ac940b58 1181 hMeanT->Reset();
1182 }// end sector loop
4d885988 1183 //return if no sector has a valid mean time
ac940b58 1184 if ( nSecMeanT == 0 ) return;
4d885988 1185
1186
7fb602b1 1187// fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1188// fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
ac940b58 1189 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1190 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1191 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1192 }
1193 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1194 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
75d8233f 1195
ac940b58 1196 fNevents++;
1197 fOldRunNumber = fRunNumber;
75d8233f 1198
ac940b58 1199 delete calIroc;
1200 delete calOroc;
1201 AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
75d8233f 1202}
1203//_____________________________________________________________________
08205ed7 1204Bool_t AliTPCCalibCE::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
1205{
1206 //
1207 // Event Processing loop - AliTPCRawStreamFast
1208 //
1209 ResetEvent();
1210 Bool_t withInput = kFALSE;
1211 while ( rawStreamFast->NextDDL() ){
4c6d06dc 1212 while ( rawStreamFast->NextChannel() ){
1213 Int_t isector = rawStreamFast->GetSector(); // current sector
1214 Int_t iRow = rawStreamFast->GetRow(); // current row
1215 Int_t iPad = rawStreamFast->GetPad(); // current pad
1216
1217 while ( rawStreamFast->NextBunch() ){
1218 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
1219 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
1220 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
1221 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
1222 Update(isector,iRow,iPad,iTimeBin+1,signal);
1223 withInput = kTRUE;
1224 }
08205ed7 1225 }
4c6d06dc 1226 }
08205ed7 1227 }
1228 if (withInput){
4c6d06dc 1229 EndEvent();
08205ed7 1230 }
1231 return withInput;
1232}
1233//_____________________________________________________________________
1234Bool_t AliTPCCalibCE::ProcessEventFast(AliRawReader *rawReader)
1235{
1236 //
1237 // Event processing loop using the fast raw stream algorithm- AliRawReader
1238 //
1239
1240 //printf("ProcessEventFast - raw reader\n");
1241
1242 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1243 if (eventHeader){
1244 fTimeStamp = eventHeader->Get("Timestamp");
1245 fRunNumber = eventHeader->Get("RunNb");
1246 }
1247 fEventId = *rawReader->GetEventId();
1248
04049c81 1249 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
08205ed7 1250 Bool_t res=ProcessEventFast(rawStreamFast);
1251 delete rawStreamFast;
1252 return res;
1253
1254}
1255//_____________________________________________________________________
7fb602b1 1256Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream)
75d8233f 1257{
1258 //
1259 // Event Processing loop - AliTPCRawStream
1260 // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
1261 //
1262
75d8233f 1263 ResetEvent();
1264
1265 Bool_t withInput = kFALSE;
1266
1267 while (rawStream->Next()) {
1268
1269 Int_t isector = rawStream->GetSector(); // current sector
1270 Int_t iRow = rawStream->GetRow(); // current row
1271 Int_t iPad = rawStream->GetPad(); // current pad
1272 Int_t iTimeBin = rawStream->GetTime(); // current time bin
1273 Float_t signal = rawStream->GetSignal(); // current ADC signal
1274
1275 Update(isector,iRow,iPad,iTimeBin,signal);
1276 withInput = kTRUE;
1277 }
1278
1279 if (withInput){
1280 EndEvent();
1281 }
1282
1283 return withInput;
1284}
1285//_____________________________________________________________________
1286Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
1287{
1288 //
1289 // Event processing loop - AliRawReader
1290 //
1291
1292
04049c81 1293 AliTPCRawStream rawStream(rawReader,(AliAltroMapping**)fMapping);
75d8233f 1294 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1295 if (eventHeader){
1296 fTimeStamp = eventHeader->Get("Timestamp");
7fb602b1 1297 fRunNumber = eventHeader->Get("RunNb");
75d8233f 1298 }
7fb602b1 1299 fEventId = *rawReader->GetEventId();
75d8233f 1300
1301
7fb602b1 1302 rawReader->Select("TPC");
75d8233f 1303
1304 return ProcessEvent(&rawStream);
1305}
1306//_____________________________________________________________________
1307Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
1308{
1309 //
1310 // Event processing loop - date event
1311 //
1312 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1313 Bool_t result=ProcessEvent(rawReader);
1314 delete rawReader;
1315 return result;
1316
1317}
1318//_____________________________________________________________________
7fb602b1 1319TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
75d8233f 1320 Int_t nbinsY, Float_t ymin, Float_t ymax,
1321 Char_t *type, Bool_t force)
1322{
1323 //
1324 // return pointer to TH2S histogram of 'type'
1325 // if force is true create a new histogram if it doesn't exist allready
1326 //
1327 if ( !force || arr->UncheckedAt(sector) )
1328 return (TH2S*)arr->UncheckedAt(sector);
1329
7fb602b1 1330 // if we are forced and histogram doesn't exist yet create it
75d8233f 1331 Char_t name[255], title[255];
1332
1333 sprintf(name,"hCalib%s%.2d",type,sector);
1334 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1335
1336 // new histogram with Q calib information. One value for each pad!
1337 TH2S* hist = new TH2S(name,title,
1338 nbinsY, ymin, ymax,
1339 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1340 hist->SetDirectory(0);
1341 arr->AddAt(hist,sector);
1342 return hist;
1343}
1344//_____________________________________________________________________
7fb602b1 1345TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
75d8233f 1346{
1347 //
1348 // return pointer to T0 histogram
1349 // if force is true create a new histogram if it doesn't exist allready
1350 //
1351 TObjArray *arr = &fHistoT0Array;
1352 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1353}
1354//_____________________________________________________________________
7fb602b1 1355TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
75d8233f 1356{
1357 //
1358 // return pointer to Q histogram
1359 // if force is true create a new histogram if it doesn't exist allready
1360 //
1361 TObjArray *arr = &fHistoQArray;
1362 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1363}
1364//_____________________________________________________________________
7fb602b1 1365TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
75d8233f 1366{
1367 //
1368 // return pointer to Q histogram
1369 // if force is true create a new histogram if it doesn't exist allready
1370 //
1371 TObjArray *arr = &fHistoRMSArray;
1372 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1373}
1374//_____________________________________________________________________
1375TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1376 Char_t *type, Bool_t force)
1377{
1378 //
1379 // return pointer to TH1S histogram
1380 // if force is true create a new histogram if it doesn't exist allready
1381 //
1382 if ( !force || arr->UncheckedAt(sector) )
1383 return (TH1S*)arr->UncheckedAt(sector);
1384
1385 // if we are forced and histogram doesn't yes exist create it
1386 Char_t name[255], title[255];
1387
1388 sprintf(name,"hCalib%s%.2d",type,sector);
1389 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1390
4d885988 1391 // new histogram with calib information. One value for each pad!
75d8233f 1392 TH1S* hist = new TH1S(name,title,
1393 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1394 hist->SetDirectory(0);
1395 arr->AddAt(hist,sector);
1396 return hist;
1397}
1398//_____________________________________________________________________
1399TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1400{
1401 //
1402 // return pointer to Q histogram
1403 // if force is true create a new histogram if it doesn't exist allready
1404 //
1405 TObjArray *arr = &fHistoTmean;
1406 return GetHisto(sector, arr, "LastTmean", force);
1407}
1408//_____________________________________________________________________
3cd27a08 1409TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
75d8233f 1410{
1411 //
1412 // return pointer to Pad Info from 'arr' for the current event and sector
1413 // if force is true create it if it doesn't exist allready
1414 //
1415 if ( !force || arr->UncheckedAt(sector) )
1416 return (TVectorF*)arr->UncheckedAt(sector);
1417
7fb602b1 1418 TVectorF *vect = new TVectorF(size);
75d8233f 1419 arr->AddAt(vect,sector);
1420 return vect;
1421}
1422//_____________________________________________________________________
7fb602b1 1423TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
75d8233f 1424{
1425 //
1426 // return pointer to Pad Times Array for the current event and sector
1427 // if force is true create it if it doesn't exist allready
1428 //
1429 TObjArray *arr = &fPadTimesArrayEvent;
7fb602b1 1430 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1431}
1432//_____________________________________________________________________
7fb602b1 1433TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
75d8233f 1434{
1435 //
1436 // return pointer to Pad Q Array for the current event and sector
1437 // if force is true create it if it doesn't exist allready
1438 // for debugging purposes only
1439 //
1440
1441 TObjArray *arr = &fPadQArrayEvent;
7fb602b1 1442 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1443}
1444//_____________________________________________________________________
7fb602b1 1445TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
75d8233f 1446{
1447 //
1448 // return pointer to Pad RMS Array for the current event and sector
1449 // if force is true create it if it doesn't exist allready
1450 // for debugging purposes only
1451 //
1452 TObjArray *arr = &fPadRMSArrayEvent;
7fb602b1 1453 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1454}
1455//_____________________________________________________________________
7fb602b1 1456TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
75d8233f 1457{
1458 //
1459 // return pointer to Pad RMS Array for the current event and sector
1460 // if force is true create it if it doesn't exist allready
1461 // for debugging purposes only
1462 //
1463 TObjArray *arr = &fPadPedestalArrayEvent;
7fb602b1 1464 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1465}
1466//_____________________________________________________________________
1467TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1468{
1469 //
1470 // return pointer to the EbyE info of the mean arrival time for 'sector'
1471 // if force is true create it if it doesn't exist allready
1472 //
1473 TObjArray *arr = &fTMeanArrayEvent;
1474 return GetVectSector(sector,arr,100,force);
75d8233f 1475}
1476//_____________________________________________________________________
7fb602b1 1477TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1478{
1479 //
1480 // return pointer to the EbyE info of the mean arrival time for 'sector'
1481 // if force is true create it if it doesn't exist allready
1482 //
1483 TObjArray *arr = &fQMeanArrayEvent;
1484 return GetVectSector(sector,arr,100,force);
1485}
1486//_____________________________________________________________________
3cd27a08 1487AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1488{
1489 //
1490 // return pointer to ROC Calibration
1491 // if force is true create a new histogram if it doesn't exist allready
1492 //
1493 if ( !force || arr->UncheckedAt(sector) )
1494 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1495
1496 // if we are forced and histogram doesn't yes exist create it
1497
1498 // new AliTPCCalROC for T0 information. One value for each pad!
1499 AliTPCCalROC *croc = new AliTPCCalROC(sector);
75d8233f 1500 arr->AddAt(croc,sector);
1501 return croc;
1502}
1503//_____________________________________________________________________
7fb602b1 1504AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
75d8233f 1505{
1506 //
ef7f7670 1507 // return pointer to Time 0 ROC Calibration
75d8233f 1508 // if force is true create a new histogram if it doesn't exist allready
1509 //
1510 TObjArray *arr = &fCalRocArrayT0;
1511 return GetCalRoc(sector, arr, force);
1512}
1513//_____________________________________________________________________
ef7f7670 1514AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1515{
1516 //
1517 // return pointer to the error of Time 0 ROC Calibration
1518 // if force is true create a new histogram if it doesn't exist allready
1519 //
1520 TObjArray *arr = &fCalRocArrayT0Err;
1521 return GetCalRoc(sector, arr, force);
1522}
1523//_____________________________________________________________________
7fb602b1 1524AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
75d8233f 1525{
1526 //
1527 // return pointer to T0 ROC Calibration
1528 // if force is true create a new histogram if it doesn't exist allready
1529 //
1530 TObjArray *arr = &fCalRocArrayQ;
1531 return GetCalRoc(sector, arr, force);
1532}
1533//_____________________________________________________________________
7fb602b1 1534AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
75d8233f 1535{
1536 //
1537 // return pointer to signal width ROC Calibration
1538 // if force is true create a new histogram if it doesn't exist allready
1539 //
1540 TObjArray *arr = &fCalRocArrayRMS;
1541 return GetCalRoc(sector, arr, force);
1542}
1543//_____________________________________________________________________
1544AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1545{
1546 //
1547 // return pointer to Outliers
1548 // if force is true create a new histogram if it doesn't exist allready
1549 //
1550 TObjArray *arr = &fCalRocArrayOutliers;
1551 return GetCalRoc(sector, arr, force);
1552}
1553//_____________________________________________________________________
3cd27a08 1554TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1555{
1556 //
1557 // return pointer to TObjArray of fit parameters
1558 // if force is true create a new histogram if it doesn't exist allready
1559 //
1560 if ( !force || arr->UncheckedAt(sector) )
1561 return (TObjArray*)arr->UncheckedAt(sector);
1562
1563 // if we are forced and array doesn't yes exist create it
1564
1565 // new TObjArray for parameters
1566 TObjArray *newArr = new TObjArray;
1567 arr->AddAt(newArr,sector);
1568 return newArr;
1569}
1570//_____________________________________________________________________
1571TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1572{
1573 //
1574 // return pointer to TObjArray of fit parameters from plane fit
1575 // if force is true create a new histogram if it doesn't exist allready
1576 //
1577 TObjArray *arr = &fParamArrayEventPol1;
1578 return GetParamArray(sector, arr, force);
1579}
1580//_____________________________________________________________________
1581TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1582{
1583 //
1584 // return pointer to TObjArray of fit parameters from parabola fit
1585 // if force is true create a new histogram if it doesn't exist allready
1586 //
1587 TObjArray *arr = &fParamArrayEventPol2;
1588 return GetParamArray(sector, arr, force);
1589}
1590//_____________________________________________________________________
7fb602b1 1591void AliTPCCalibCE::ResetEvent()
75d8233f 1592{
1593 //
1594 // Reset global counters -- Should be called before each event is processed
1595 //
1596 fLastSector=-1;
1597 fCurrentSector=-1;
1598 fCurrentRow=-1;
1599 fCurrentChannel=-1;
1600
1601 ResetPad();
1602
1603 fPadTimesArrayEvent.Delete();
1604 fPadQArrayEvent.Delete();
1605 fPadRMSArrayEvent.Delete();
1606 fPadPedestalArrayEvent.Delete();
1607
7fb602b1 1608 for ( Int_t i=0; i<72; ++i ){
bf57d87d 1609 fVTime0Offset.GetMatrixArray()[i]=0;
1610 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1611 fVMeanQ.GetMatrixArray()[i]=0;
1612 fVMeanQCounter.GetMatrixArray()[i]=0;
75d8233f 1613 }
1614}
1615//_____________________________________________________________________
7fb602b1 1616void AliTPCCalibCE::ResetPad()
75d8233f 1617{
1618 //
1619 // Reset pad infos -- Should be called after a pad has been processed
1620 //
7fb602b1 1621 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
bf57d87d 1622 fPadSignal.GetMatrixArray()[i] = 0;
75d8233f 1623 fMaxTimeBin = -1;
1624 fMaxPadSignal = -1;
1625 fPadPedestal = -1;
1626 fPadNoise = -1;
1627}
1628//_____________________________________________________________________
7fb602b1 1629void AliTPCCalibCE::Merge(AliTPCCalibCE *ce)
1630{
1631 //
1632 // Merge ce to the current AliTPCCalibCE
1633 //
1634
1635 //merge histograms
1636 for (Int_t iSec=0; iSec<72; ++iSec){
1637 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1638 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1639 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1640
1641
1642 if ( hRefQmerge ){
1643 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1644 TH2S *hRefQ = GetHistoQ(iSec);
1645 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1646 else {
1647 TH2S *hist = new TH2S(*hRefQmerge);
1648 hist->SetDirectory(0);
1649 fHistoQArray.AddAt(hist, iSec);
1650 }
1651 hRefQmerge->SetDirectory(dir);
1652 }
1653 if ( hRefT0merge ){
1654 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1655 TH2S *hRefT0 = GetHistoT0(iSec);
1656 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1657 else {
1658 TH2S *hist = new TH2S(*hRefT0merge);
1659 hist->SetDirectory(0);
1660 fHistoT0Array.AddAt(hist, iSec);
1661 }
1662 hRefT0merge->SetDirectory(dir);
1663 }
1664 if ( hRefRMSmerge ){
1665 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1666 TH2S *hRefRMS = GetHistoRMS(iSec);
1667 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1668 else {
1669 TH2S *hist = new TH2S(*hRefRMSmerge);
1670 hist->SetDirectory(0);
1671 fHistoRMSArray.AddAt(hist, iSec);
1672 }
1673 hRefRMSmerge->SetDirectory(dir);
1674 }
1675
1676 }
1677
1678 // merge time information
1679
1680
1681 Int_t nCEevents = ce->GetNeventsProcessed();
1682 for (Int_t iSec=0; iSec<72; ++iSec){
1683 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1684 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1685 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1686 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1687
1688 TObjArray *arrPol1 = 0x0;
1689 TObjArray *arrPol2 = 0x0;
1690 TVectorF *vMeanTime = 0x0;
1691 TVectorF *vMeanQ = 0x0;
1692
1693 //resize arrays
1694 if ( arrPol1CE && arrPol2CE ){
1695 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1696 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1697 arrPol1->Expand(fNevents+nCEevents);
1698 arrPol2->Expand(fNevents+nCEevents);
1699 }
1700 if ( vMeanTimeCE && vMeanQCE ){
d8faeea7 1701 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1702 vMeanQ = GetQMeanEvents(iSec,kTRUE);
7fb602b1 1703 vMeanTime->ResizeTo(fNevents+nCEevents);
1704 vMeanQ->ResizeTo(fNevents+nCEevents);
1705 }
1706
1707
1708 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1709 if ( arrPol1CE && arrPol2CE ){
1710 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1711 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1712 if ( paramPol1 && paramPol2 ){
1713 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1714 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1715 }
1716 }
1717 if ( vMeanTimeCE && vMeanQCE ){
1718 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1719 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1720 }
1721 }
1722 }
1723
1724
1725
1726 TVectorD* eventTimes = ce->GetEventTimes();
1727 TVectorD* eventIds = ce->GetEventIds();
1728 fVEventTime.ResizeTo(fNevents+nCEevents);
1729 fVEventNumber.ResizeTo(fNevents+nCEevents);
1730
1731 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1732 Double_t evTime = eventTimes->GetMatrixArray()[iEvent];
1733 Double_t evId = eventIds->GetMatrixArray()[iEvent];
1734 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1735 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1736 }
1737 fNevents+=nCEevents; //increase event counter
1738
1739}
1740//_____________________________________________________________________
1741TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
75d8233f 1742{
ac940b58 1743 //
1744 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1745 // xVariable: 0-event time, 1-event id, 2-internal event counter
1746 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1747 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1748 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1749 // not used for mean time and mean Q )
1750 // for an example see class description at the beginning
1751 //
75d8233f 1752
ac940b58 1753 Double_t *x = new Double_t[fNevents];
1754 Double_t *y = new Double_t[fNevents];
75d8233f 1755
ac940b58 1756 TVectorD *xVar = 0x0;
1757 TObjArray *aType = 0x0;
1758 Int_t npoints=0;
75d8233f 1759
1760 // sanity checks
ac940b58 1761 if ( !GetHistoT0(sector) ) return 0x0; //Sector has not been filled
1762 if ( (sector<0) || (sector>71) ) return 0x0;
1763 if ( (xVariable<0) || (xVariable>2) ) return 0x0;
1764 if ( (fitType<0) || (fitType>3) ) return 0x0;
1765 if ( !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1766
1767 if ( fitType==0 ){
1768 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1769 aType = &fParamArrayEventPol1;
1770 if ( aType->At(sector)==0x0 ) return 0x0;
1771 }
1772 else if ( fitType==1 ){
1773 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1774 aType = &fParamArrayEventPol2;
1775 if ( aType->At(sector)==0x0 ) return 0x0;
1776 }
75d8233f 1777
1778
ac940b58 1779 if ( xVariable == 0 ) xVar = &fVEventTime;
1780 if ( xVariable == 1 ) xVar = &fVEventNumber;
1781 if ( xVariable == 2 ) {
1782 xVar = new TVectorD(fNevents);
1783 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1784 }
75d8233f 1785
ac940b58 1786 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1787 if ( fitType<2 ){
1788 TObjArray *events = (TObjArray*)(aType->At(sector));
1789 if ( events->GetSize()<=ievent ) break;
1790 TVectorD *v = (TVectorD*)(events->At(ievent));
1791 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1792 } else if (fitType == 2) {
1793 Double_t xValue=(*xVar)[ievent];
1794 Double_t yValue=(*GetTMeanEvents(sector))[ievent];
1795 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1796 }else if (fitType == 3) {
1797 Double_t xValue=(*xVar)[ievent];
1798 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1799 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
75d8233f 1800 }
ac940b58 1801 }
75d8233f 1802
ac940b58 1803 TGraph *gr = new TGraph(npoints);
75d8233f 1804 //sort xVariable increasing
ac940b58 1805 Int_t *sortIndex = new Int_t[npoints];
1806 TMath::Sort(npoints,x,sortIndex);
1807 for (Int_t i=0;i<npoints;++i){
1808 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1809 }
75d8233f 1810
1811
ac940b58 1812 if ( xVariable == 2 ) delete xVar;
1813 delete x;
1814 delete y;
1815 delete sortIndex;
1816 return gr;
75d8233f 1817}
1818//_____________________________________________________________________
1819void AliTPCCalibCE::Analyse()
1820{
1821 //
1822 // Calculate calibration constants
1823 //
1824
1825 TVectorD paramQ(3);
1826 TVectorD paramT0(3);
1827 TVectorD paramRMS(3);
1828 TMatrixD dummy(3,3);
1829
ef7f7670 1830 Float_t channelCounter=0;
1831 fMeanT0rms=0;
1832 fMeanQrms=0;
1833 fMeanRMSrms=0;
1834
7fb602b1 1835 for (Int_t iSec=0; iSec<72; ++iSec){
75d8233f 1836 TH2S *hT0 = GetHistoT0(iSec);
1837 if (!hT0 ) continue;
1838
ef7f7670 1839 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1840 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1841 AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1842 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1843 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
75d8233f 1844
1845 TH2S *hQ = GetHistoQ(iSec);
1846 TH2S *hRMS = GetHistoRMS(iSec);
1847
3cd27a08 1848 Short_t *arrayhQ = hQ->GetArray();
1849 Short_t *arrayhT0 = hT0->GetArray();
1850 Short_t *arrayhRMS = hRMS->GetArray();
75d8233f 1851
1852 UInt_t nChannels = fROC->GetNChannels(iSec);
1853
1854 //debug
1855 Int_t row=0;
1856 Int_t pad=0;
1857 Int_t padc=0;
1858 //! debug
1859
7fb602b1 1860 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
75d8233f 1861
1862
1863 Float_t cogTime0 = -1000;
1864 Float_t cogQ = -1000;
1865 Float_t cogRMS = -1000;
ef7f7670 1866 Float_t cogOut = 0;
1867 Float_t rms = 0;
1868 Float_t rmsT0 = 0;
75d8233f 1869
1870
1871 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1872 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1873 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1874
ef7f7670 1875 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
1876 fMeanQrms+=rms;
1877 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
1878 fMeanT0rms+=rmsT0;
1879 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
1880 fMeanRMSrms+=rms;
1881 channelCounter++;
75d8233f 1882
1883 /*
7fb602b1 1884 //outlier specifications
75d8233f 1885 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1886 cogOut = 1;
1887 cogTime0 = 0;
1888 cogQ = 0;
1889 cogRMS = 0;
1890 }
1891*/
1892 rocQ->SetValue(iChannel, cogQ*cogQ);
1893 rocT0->SetValue(iChannel, cogTime0);
ef7f7670 1894 rocT0Err->SetValue(iChannel, rmsT0);
75d8233f 1895 rocRMS->SetValue(iChannel, cogRMS);
1896 rocOut->SetValue(iChannel, cogOut);
1897
1898
1899 //debug
1900 if ( fDebugLevel > 0 ){
1901 if ( !fDebugStreamer ) {
1902 //debug stream
1903 TDirectory *backup = gDirectory;
1904 fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1905 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1906 }
1907
1908 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1909 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1910 padc = pad-(fROC->GetNPads(iSec,row)/2);
1911
1912 (*fDebugStreamer) << "DataEnd" <<
1913 "Sector=" << iSec <<
1914 "Pad=" << pad <<
1915 "PadC=" << padc <<
1916 "Row=" << row <<
1917 "PadSec=" << iChannel <<
1918 "Q=" << cogQ <<
1919 "T0=" << cogTime0 <<
1920 "RMS=" << cogRMS <<
1921 "\n";
1922 }
1923 //! debug
1924
1925 }
1926
1927 }
ef7f7670 1928 if ( channelCounter>0 ){
1929 fMeanT0rms/=channelCounter;
1930 fMeanQrms/=channelCounter;
1931 fMeanRMSrms/=channelCounter;
1932 }
7fb602b1 1933 if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
75d8233f 1934// delete fDebugStreamer;
bf57d87d 1935// fDebugStreamer = 0x0;
75d8233f 1936}
1937//_____________________________________________________________________
1938void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1939{
1940 //
1941 // Write class to file
1942 //
1943
1944 TString sDir(dir);
1945 TString option;
1946
1947 if ( append )
1948 option = "update";
1949 else
1950 option = "recreate";
1951
1952 TDirectory *backup = gDirectory;
1953 TFile f(filename,option.Data());
1954 f.cd();
1955 if ( !sDir.IsNull() ){
1956 f.mkdir(sDir.Data());
1957 f.cd(sDir);
1958 }
1959 this->Write();
1960 f.Close();
1961
1962 if ( backup ) backup->cd();
1963}