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