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