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