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