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