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