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