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