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