]> 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
Possibility to getthe data from zip file
[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
75d8233f 852 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
7fb602b1 853 for ( Int_t iSec = 0; iSec<72; ++iSec ){
bf57d87d 854 //find median and then calculate the mean around it
7fb602b1 855 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
bf57d87d 856 if ( !hMeanT ) continue;
7fb602b1 857
75d8233f 858 Double_t entries = hMeanT->GetEntries();
859 Double_t sum = 0;
860 Short_t *arr = hMeanT->GetArray()+1;
861 Int_t ibin=0;
7fb602b1 862 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
75d8233f 863 sum+=arr[ibin];
864 if ( sum>=(entries/2) ) break;
865 }
866 Int_t delta = 4;
867 Int_t firstBin = fFirstTimeBin+ibin-delta;
868 Int_t lastBin = fFirstTimeBin+ibin+delta;
869 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
870 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
bf57d87d 871 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
7fb602b1 872
873 // check boundaries for ebye info of mean time
874 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
875 Int_t vSize=vMeanTime->GetNrows();
876 if ( vSize < fNevents+1 )
877 vMeanTime->ResizeTo(vSize+100);
878
879 vMeanTime->GetMatrixArray()[fNevents]=median;
bf57d87d 880 // end find median
881
882 TVectorF *vTimes = GetPadTimesEvent(iSec);
7fb602b1 883 if ( !vTimes ) continue; //continue if no time information for this sector is available
884
885
bf57d87d 886 AliTPCCalROC calIrocOutliers(0);
887 AliTPCCalROC calOrocOutliers(36);
888
889 // calculate mean Q of the sector
890 Float_t meanQ = 0;
891 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
7fb602b1 892 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
893 if ( vSize < fNevents+1 ) // vSize is the same as for vMeanTime!
894 vMeanQ->ResizeTo(vSize+100);
75d8233f 895
7fb602b1 896 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
897
898 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
3cd27a08 899 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
75d8233f 900
901 //set values for temporary roc calibration class
bf57d87d 902 if ( iSec < 36 ) {
3cd27a08 903 calIroc->SetValue(iChannel, time);
904 if ( time == 0 ) calIrocOutliers.SetValue(iChannel,1);
bf57d87d 905
906 } else {
3cd27a08 907 calOroc->SetValue(iChannel, time);
908 if ( time == 0 ) calOrocOutliers.SetValue(iChannel,1);
bf57d87d 909 }
75d8233f 910
911 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
3cd27a08 912 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
75d8233f 913
914
915
916 //------------------------------- Debug start ------------------------------
917 if ( fDebugLevel>0 ){
918 if ( !fDebugStreamer ) {
919 //debug stream
920 TDirectory *backup = gDirectory;
921 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
922 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
923 }
924
925 Int_t row=0;
926 Int_t pad=0;
927 Int_t padc=0;
928
3cd27a08 929 Float_t q = (*GetPadQEvent(iSec))[iChannel];
930 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
75d8233f 931
932 UInt_t channel=iChannel;
933 Int_t sector=iSec;
934
935 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
936 pad = channel-fROC->GetRowIndexes(sector)[row];
937 padc = pad-(fROC->GetNPads(sector,row)/2);
938
939// TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
940// Form("hSignalD%d.%d.%d",sector,row,pad),
941// fLastTimeBin-fFirstTimeBin,
942// fFirstTimeBin,fLastTimeBin);
943// h1->SetDirectory(0);
944//
7fb602b1 945// for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
75d8233f 946// h1->Fill(i,fPadSignal(i));
947
3cd27a08 948 Double_t t0Sec = 0;
bf57d87d 949 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
3cd27a08 950 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
951 Double_t t0Side = time0Side[(iSec/18)%2];
75d8233f 952 (*fDebugStreamer) << "DataPad" <<
953 "Event=" << fNevents <<
7fb602b1 954 "RunNumber=" << fRunNumber <<
bf57d87d 955 "TimeStamp=" << fTimeStamp <<
75d8233f 956 "Sector="<< sector <<
957 "Row=" << row<<
958 "Pad=" << pad <<
959 "PadC=" << padc <<
960 "PadSec="<< channel <<
3cd27a08 961 "Time0Sec=" << t0Sec <<
962 "Time0Side=" << t0Side <<
963 "Time=" << time <<
964 "RMS=" << rms <<
965 "Sum=" << q <<
bf57d87d 966 "MeanQ=" << meanQ <<
967 // "hist.=" << h1 <<
75d8233f 968 "\n";
969
bf57d87d 970 // delete h1;
971
75d8233f 972 }
973 //----------------------------- Debug end ------------------------------
974 }// end channel loop
975 hMeanT->Reset();
976
977 TVectorD paramPol1(3);
978 TVectorD paramPol2(6);
979 TMatrixD matPol1(3,3);
980 TMatrixD matPol2(6,6);
bf57d87d 981 Float_t chi2Pol1=0;
982 Float_t chi2Pol2=0;
75d8233f 983
984 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
985 if ( iSec < 36 ){
7fb602b1 986 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
987 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
75d8233f 988 } else {
7fb602b1 989 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
990 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
75d8233f 991 }
992
993 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
994 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
995 }
bf57d87d 996// printf("events: %d -- size: %d\n",fNevents,GetParamArrayPol1(iSec)->GetSize());
75d8233f 997
998 //------------------------------- Debug start ------------------------------
999 if ( fDebugLevel>0 ){
bf57d87d 1000 if ( !fDebugStreamer ) {
1001 //debug stream
1002 TDirectory *backup = gDirectory;
1003 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
1004 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1005 }
75d8233f 1006 (*fDebugStreamer) << "DataRoc" <<
3cd27a08 1007// "Event=" << fEvent <<
7fb602b1 1008 "RunNumber=" << fRunNumber <<
75d8233f 1009 "TimeStamp=" << fTimeStamp <<
1010 "Sector="<< iSec <<
1011 "hMeanT.=" << hMeanT <<
1012 "median=" << median <<
1013 "paramPol1.=" << &paramPol1 <<
1014 "paramPol2.=" << &paramPol2 <<
1015 "matPol1.=" << &matPol1 <<
1016 "matPol2.=" << &matPol2 <<
1017 "chi2Pol1=" << chi2Pol1 <<
1018 "chi2Pol2=" << chi2Pol2 <<
1019 "\n";
1020 }
1021 //------------------------------- Debug end ------------------------------
1022 }// end sector loop
1023
7fb602b1 1024// fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1025// fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1026 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1027 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1028 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
75d8233f 1029 }
7fb602b1 1030 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1031 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
75d8233f 1032
1033 fNevents++;
1034 fOldRunNumber = fRunNumber;
1035
7fb602b1 1036 delete calIroc;
1037 delete calOroc;
75d8233f 1038}
1039//_____________________________________________________________________
08205ed7 1040Bool_t AliTPCCalibCE::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
1041{
1042 //
1043 // Event Processing loop - AliTPCRawStreamFast
1044 //
1045 ResetEvent();
1046 Bool_t withInput = kFALSE;
1047 while ( rawStreamFast->NextDDL() ){
1048 while ( rawStreamFast->NextChannel() ){
1049 Int_t isector = rawStreamFast->GetSector(); // current sector
1050 Int_t iRow = rawStreamFast->GetRow(); // current row
1051 Int_t iPad = rawStreamFast->GetPad(); // current pad
1052 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
1053 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
1054
1055 while ( rawStreamFast->NextBunch() ){
1056 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
1057 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
1058 Update(isector,iRow,iPad,iTimeBin+1,signal);
1059 withInput = kTRUE;
1060 }
1061 }
1062 }
1063 }
1064 if (withInput){
1065 EndEvent();
1066 }
1067 return withInput;
1068}
1069//_____________________________________________________________________
1070Bool_t AliTPCCalibCE::ProcessEventFast(AliRawReader *rawReader)
1071{
1072 //
1073 // Event processing loop using the fast raw stream algorithm- AliRawReader
1074 //
1075
1076 //printf("ProcessEventFast - raw reader\n");
1077
1078 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1079 if (eventHeader){
1080 fTimeStamp = eventHeader->Get("Timestamp");
1081 fRunNumber = eventHeader->Get("RunNb");
1082 }
1083 fEventId = *rawReader->GetEventId();
1084
04049c81 1085 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader, (AliAltroMapping**)fMapping);
08205ed7 1086 Bool_t res=ProcessEventFast(rawStreamFast);
1087 delete rawStreamFast;
1088 return res;
1089
1090}
1091//_____________________________________________________________________
7fb602b1 1092Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream)
75d8233f 1093{
1094 //
1095 // Event Processing loop - AliTPCRawStream
1096 // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
1097 //
1098
1099 rawStream->SetOldRCUFormat(fOldRCUformat);
1100
1101 ResetEvent();
1102
1103 Bool_t withInput = kFALSE;
1104
1105 while (rawStream->Next()) {
1106
1107 Int_t isector = rawStream->GetSector(); // current sector
1108 Int_t iRow = rawStream->GetRow(); // current row
1109 Int_t iPad = rawStream->GetPad(); // current pad
1110 Int_t iTimeBin = rawStream->GetTime(); // current time bin
1111 Float_t signal = rawStream->GetSignal(); // current ADC signal
1112
1113 Update(isector,iRow,iPad,iTimeBin,signal);
1114 withInput = kTRUE;
1115 }
1116
1117 if (withInput){
1118 EndEvent();
1119 }
1120
1121 return withInput;
1122}
1123//_____________________________________________________________________
1124Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
1125{
1126 //
1127 // Event processing loop - AliRawReader
1128 //
1129
1130
04049c81 1131 AliTPCRawStream rawStream(rawReader,(AliAltroMapping**)fMapping);
75d8233f 1132 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
1133 if (eventHeader){
1134 fTimeStamp = eventHeader->Get("Timestamp");
7fb602b1 1135 fRunNumber = eventHeader->Get("RunNb");
75d8233f 1136 }
7fb602b1 1137 fEventId = *rawReader->GetEventId();
75d8233f 1138
1139
7fb602b1 1140 rawReader->Select("TPC");
75d8233f 1141
1142 return ProcessEvent(&rawStream);
1143}
1144//_____________________________________________________________________
1145Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
1146{
1147 //
1148 // Event processing loop - date event
1149 //
1150 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
1151 Bool_t result=ProcessEvent(rawReader);
1152 delete rawReader;
1153 return result;
1154
1155}
1156//_____________________________________________________________________
7fb602b1 1157TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
75d8233f 1158 Int_t nbinsY, Float_t ymin, Float_t ymax,
1159 Char_t *type, Bool_t force)
1160{
1161 //
1162 // return pointer to TH2S histogram of 'type'
1163 // if force is true create a new histogram if it doesn't exist allready
1164 //
1165 if ( !force || arr->UncheckedAt(sector) )
1166 return (TH2S*)arr->UncheckedAt(sector);
1167
7fb602b1 1168 // if we are forced and histogram doesn't exist yet create it
75d8233f 1169 Char_t name[255], title[255];
1170
1171 sprintf(name,"hCalib%s%.2d",type,sector);
1172 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1173
1174 // new histogram with Q calib information. One value for each pad!
1175 TH2S* hist = new TH2S(name,title,
1176 nbinsY, ymin, ymax,
1177 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1178 hist->SetDirectory(0);
1179 arr->AddAt(hist,sector);
1180 return hist;
1181}
1182//_____________________________________________________________________
7fb602b1 1183TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
75d8233f 1184{
1185 //
1186 // return pointer to T0 histogram
1187 // if force is true create a new histogram if it doesn't exist allready
1188 //
1189 TObjArray *arr = &fHistoT0Array;
1190 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1191}
1192//_____________________________________________________________________
7fb602b1 1193TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
75d8233f 1194{
1195 //
1196 // return pointer to Q histogram
1197 // if force is true create a new histogram if it doesn't exist allready
1198 //
1199 TObjArray *arr = &fHistoQArray;
1200 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1201}
1202//_____________________________________________________________________
7fb602b1 1203TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
75d8233f 1204{
1205 //
1206 // return pointer to Q histogram
1207 // if force is true create a new histogram if it doesn't exist allready
1208 //
1209 TObjArray *arr = &fHistoRMSArray;
1210 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1211}
1212//_____________________________________________________________________
1213TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1214 Char_t *type, Bool_t force)
1215{
1216 //
1217 // return pointer to TH1S histogram
1218 // if force is true create a new histogram if it doesn't exist allready
1219 //
1220 if ( !force || arr->UncheckedAt(sector) )
1221 return (TH1S*)arr->UncheckedAt(sector);
1222
1223 // if we are forced and histogram doesn't yes exist create it
1224 Char_t name[255], title[255];
1225
1226 sprintf(name,"hCalib%s%.2d",type,sector);
1227 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
1228
1229 // new histogram with Q calib information. One value for each pad!
1230 TH1S* hist = new TH1S(name,title,
1231 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1232 hist->SetDirectory(0);
1233 arr->AddAt(hist,sector);
1234 return hist;
1235}
1236//_____________________________________________________________________
1237TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1238{
1239 //
1240 // return pointer to Q histogram
1241 // if force is true create a new histogram if it doesn't exist allready
1242 //
1243 TObjArray *arr = &fHistoTmean;
1244 return GetHisto(sector, arr, "LastTmean", force);
1245}
1246//_____________________________________________________________________
3cd27a08 1247TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
75d8233f 1248{
1249 //
1250 // return pointer to Pad Info from 'arr' for the current event and sector
1251 // if force is true create it if it doesn't exist allready
1252 //
1253 if ( !force || arr->UncheckedAt(sector) )
1254 return (TVectorF*)arr->UncheckedAt(sector);
1255
7fb602b1 1256 TVectorF *vect = new TVectorF(size);
75d8233f 1257 arr->AddAt(vect,sector);
1258 return vect;
1259}
1260//_____________________________________________________________________
7fb602b1 1261TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
75d8233f 1262{
1263 //
1264 // return pointer to Pad Times Array for the current event and sector
1265 // if force is true create it if it doesn't exist allready
1266 //
1267 TObjArray *arr = &fPadTimesArrayEvent;
7fb602b1 1268 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1269}
1270//_____________________________________________________________________
7fb602b1 1271TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
75d8233f 1272{
1273 //
1274 // return pointer to Pad Q Array for the current event and sector
1275 // if force is true create it if it doesn't exist allready
1276 // for debugging purposes only
1277 //
1278
1279 TObjArray *arr = &fPadQArrayEvent;
7fb602b1 1280 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1281}
1282//_____________________________________________________________________
7fb602b1 1283TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
75d8233f 1284{
1285 //
1286 // return pointer to Pad RMS Array for the current event and sector
1287 // if force is true create it if it doesn't exist allready
1288 // for debugging purposes only
1289 //
1290 TObjArray *arr = &fPadRMSArrayEvent;
7fb602b1 1291 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1292}
1293//_____________________________________________________________________
7fb602b1 1294TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
75d8233f 1295{
1296 //
1297 // return pointer to Pad RMS Array for the current event and sector
1298 // if force is true create it if it doesn't exist allready
1299 // for debugging purposes only
1300 //
1301 TObjArray *arr = &fPadPedestalArrayEvent;
7fb602b1 1302 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1303}
1304//_____________________________________________________________________
1305TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1306{
1307 //
1308 // return pointer to the EbyE info of the mean arrival time for 'sector'
1309 // if force is true create it if it doesn't exist allready
1310 //
1311 TObjArray *arr = &fTMeanArrayEvent;
1312 return GetVectSector(sector,arr,100,force);
75d8233f 1313}
1314//_____________________________________________________________________
7fb602b1 1315TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1316{
1317 //
1318 // return pointer to the EbyE info of the mean arrival time for 'sector'
1319 // if force is true create it if it doesn't exist allready
1320 //
1321 TObjArray *arr = &fQMeanArrayEvent;
1322 return GetVectSector(sector,arr,100,force);
1323}
1324//_____________________________________________________________________
3cd27a08 1325AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1326{
1327 //
1328 // return pointer to ROC Calibration
1329 // if force is true create a new histogram if it doesn't exist allready
1330 //
1331 if ( !force || arr->UncheckedAt(sector) )
1332 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1333
1334 // if we are forced and histogram doesn't yes exist create it
1335
1336 // new AliTPCCalROC for T0 information. One value for each pad!
1337 AliTPCCalROC *croc = new AliTPCCalROC(sector);
75d8233f 1338 arr->AddAt(croc,sector);
1339 return croc;
1340}
1341//_____________________________________________________________________
7fb602b1 1342AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
75d8233f 1343{
1344 //
1345 // return pointer to Carge ROC Calibration
1346 // if force is true create a new histogram if it doesn't exist allready
1347 //
1348 TObjArray *arr = &fCalRocArrayT0;
1349 return GetCalRoc(sector, arr, force);
1350}
1351//_____________________________________________________________________
7fb602b1 1352AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
75d8233f 1353{
1354 //
1355 // return pointer to T0 ROC Calibration
1356 // if force is true create a new histogram if it doesn't exist allready
1357 //
1358 TObjArray *arr = &fCalRocArrayQ;
1359 return GetCalRoc(sector, arr, force);
1360}
1361//_____________________________________________________________________
7fb602b1 1362AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
75d8233f 1363{
1364 //
1365 // return pointer to signal width ROC Calibration
1366 // if force is true create a new histogram if it doesn't exist allready
1367 //
1368 TObjArray *arr = &fCalRocArrayRMS;
1369 return GetCalRoc(sector, arr, force);
1370}
1371//_____________________________________________________________________
1372AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1373{
1374 //
1375 // return pointer to Outliers
1376 // if force is true create a new histogram if it doesn't exist allready
1377 //
1378 TObjArray *arr = &fCalRocArrayOutliers;
1379 return GetCalRoc(sector, arr, force);
1380}
1381//_____________________________________________________________________
3cd27a08 1382TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1383{
1384 //
1385 // return pointer to TObjArray of fit parameters
1386 // if force is true create a new histogram if it doesn't exist allready
1387 //
1388 if ( !force || arr->UncheckedAt(sector) )
1389 return (TObjArray*)arr->UncheckedAt(sector);
1390
1391 // if we are forced and array doesn't yes exist create it
1392
1393 // new TObjArray for parameters
1394 TObjArray *newArr = new TObjArray;
1395 arr->AddAt(newArr,sector);
1396 return newArr;
1397}
1398//_____________________________________________________________________
1399TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1400{
1401 //
1402 // return pointer to TObjArray of fit parameters from plane fit
1403 // if force is true create a new histogram if it doesn't exist allready
1404 //
1405 TObjArray *arr = &fParamArrayEventPol1;
1406 return GetParamArray(sector, arr, force);
1407}
1408//_____________________________________________________________________
1409TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1410{
1411 //
1412 // return pointer to TObjArray of fit parameters from parabola fit
1413 // if force is true create a new histogram if it doesn't exist allready
1414 //
1415 TObjArray *arr = &fParamArrayEventPol2;
1416 return GetParamArray(sector, arr, force);
1417}
1418//_____________________________________________________________________
7fb602b1 1419void AliTPCCalibCE::ResetEvent()
75d8233f 1420{
1421 //
1422 // Reset global counters -- Should be called before each event is processed
1423 //
1424 fLastSector=-1;
1425 fCurrentSector=-1;
1426 fCurrentRow=-1;
1427 fCurrentChannel=-1;
1428
1429 ResetPad();
1430
1431 fPadTimesArrayEvent.Delete();
1432 fPadQArrayEvent.Delete();
1433 fPadRMSArrayEvent.Delete();
1434 fPadPedestalArrayEvent.Delete();
1435
7fb602b1 1436 for ( Int_t i=0; i<72; ++i ){
bf57d87d 1437 fVTime0Offset.GetMatrixArray()[i]=0;
1438 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1439 fVMeanQ.GetMatrixArray()[i]=0;
1440 fVMeanQCounter.GetMatrixArray()[i]=0;
75d8233f 1441 }
1442}
1443//_____________________________________________________________________
7fb602b1 1444void AliTPCCalibCE::ResetPad()
75d8233f 1445{
1446 //
1447 // Reset pad infos -- Should be called after a pad has been processed
1448 //
7fb602b1 1449 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
bf57d87d 1450 fPadSignal.GetMatrixArray()[i] = 0;
75d8233f 1451 fMaxTimeBin = -1;
1452 fMaxPadSignal = -1;
1453 fPadPedestal = -1;
1454 fPadNoise = -1;
1455}
1456//_____________________________________________________________________
7fb602b1 1457void AliTPCCalibCE::Merge(AliTPCCalibCE *ce)
1458{
1459 //
1460 // Merge ce to the current AliTPCCalibCE
1461 //
1462
1463 //merge histograms
1464 for (Int_t iSec=0; iSec<72; ++iSec){
1465 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1466 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1467 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1468
1469
1470 if ( hRefQmerge ){
1471 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1472 TH2S *hRefQ = GetHistoQ(iSec);
1473 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1474 else {
1475 TH2S *hist = new TH2S(*hRefQmerge);
1476 hist->SetDirectory(0);
1477 fHistoQArray.AddAt(hist, iSec);
1478 }
1479 hRefQmerge->SetDirectory(dir);
1480 }
1481 if ( hRefT0merge ){
1482 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1483 TH2S *hRefT0 = GetHistoT0(iSec);
1484 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1485 else {
1486 TH2S *hist = new TH2S(*hRefT0merge);
1487 hist->SetDirectory(0);
1488 fHistoT0Array.AddAt(hist, iSec);
1489 }
1490 hRefT0merge->SetDirectory(dir);
1491 }
1492 if ( hRefRMSmerge ){
1493 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1494 TH2S *hRefRMS = GetHistoRMS(iSec);
1495 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1496 else {
1497 TH2S *hist = new TH2S(*hRefRMSmerge);
1498 hist->SetDirectory(0);
1499 fHistoRMSArray.AddAt(hist, iSec);
1500 }
1501 hRefRMSmerge->SetDirectory(dir);
1502 }
1503
1504 }
1505
1506 // merge time information
1507
1508
1509 Int_t nCEevents = ce->GetNeventsProcessed();
1510 for (Int_t iSec=0; iSec<72; ++iSec){
1511 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1512 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1513 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1514 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1515
1516 TObjArray *arrPol1 = 0x0;
1517 TObjArray *arrPol2 = 0x0;
1518 TVectorF *vMeanTime = 0x0;
1519 TVectorF *vMeanQ = 0x0;
1520
1521 //resize arrays
1522 if ( arrPol1CE && arrPol2CE ){
1523 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1524 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1525 arrPol1->Expand(fNevents+nCEevents);
1526 arrPol2->Expand(fNevents+nCEevents);
1527 }
1528 if ( vMeanTimeCE && vMeanQCE ){
d8faeea7 1529 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1530 vMeanQ = GetQMeanEvents(iSec,kTRUE);
7fb602b1 1531 vMeanTime->ResizeTo(fNevents+nCEevents);
1532 vMeanQ->ResizeTo(fNevents+nCEevents);
1533 }
1534
1535
1536 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1537 if ( arrPol1CE && arrPol2CE ){
1538 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1539 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1540 if ( paramPol1 && paramPol2 ){
1541 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1542 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1543 }
1544 }
1545 if ( vMeanTimeCE && vMeanQCE ){
1546 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1547 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1548 }
1549 }
1550 }
1551
1552
1553
1554 TVectorD* eventTimes = ce->GetEventTimes();
1555 TVectorD* eventIds = ce->GetEventIds();
1556 fVEventTime.ResizeTo(fNevents+nCEevents);
1557 fVEventNumber.ResizeTo(fNevents+nCEevents);
1558
1559 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1560 Double_t evTime = eventTimes->GetMatrixArray()[iEvent];
1561 Double_t evId = eventIds->GetMatrixArray()[iEvent];
1562 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1563 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1564 }
1565 fNevents+=nCEevents; //increase event counter
1566
1567}
1568//_____________________________________________________________________
1569TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
75d8233f 1570{
1571 //
7fb602b1 1572 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1573 // xVariable: 0-event time, 1-event id, 2-internal event counter
1574 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1575 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1576 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1577 // not used for mean time and mean Q )
1578 // for an example see class description at the beginning
75d8233f 1579 //
1580
1581 Double_t *x = new Double_t[fNevents];
1582 Double_t *y = new Double_t[fNevents];
1583
1584 TVectorD *xVar = 0x0;
1585 TObjArray *aType = 0x0;
1586 Int_t npoints=0;
1587
1588 // sanity checks
1589 if ( (sector<0) || (sector>71) ) return 0x0;
1590 if ( (xVariable<0) || (xVariable>2) ) return 0x0;
bf57d87d 1591 if ( (fitType<0) || (fitType>3) ) return 0x0;
75d8233f 1592 if ( fitType==0 ){
1593 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1594 aType = &fParamArrayEventPol1;
1595 if ( aType->At(sector)==0x0 ) return 0x0;
1596 }
1597 else if ( fitType==1 ){
1598 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1599 aType = &fParamArrayEventPol2;
1600 if ( aType->At(sector)==0x0 ) return 0x0;
1601 }
1602
1603
1604 if ( xVariable == 0 ) xVar = &fVEventTime;
1605 if ( xVariable == 1 ) xVar = &fVEventNumber;
1606 if ( xVariable == 2 ) {
1607 xVar = new TVectorD(fNevents);
7fb602b1 1608 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
75d8233f 1609 }
1610
7fb602b1 1611 for (Int_t ievent =0; ievent<fNevents; ++ievent){
75d8233f 1612 if ( fitType<2 ){
1613 TObjArray *events = (TObjArray*)(aType->At(sector));
1614 if ( events->GetSize()<=ievent ) break;
1615 TVectorD *v = (TVectorD*)(events->At(ievent));
7fb602b1 1616 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
75d8233f 1617 } else if (fitType == 2) {
7fb602b1 1618 Double_t xValue=(*xVar)[ievent];
1619 Double_t yValue=(*GetTMeanEvents(sector))[ievent];
1620 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
bf57d87d 1621 }else if (fitType == 3) {
7fb602b1 1622 Double_t xValue=(*xVar)[ievent];
1623 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1624 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
75d8233f 1625 }
1626 }
1627
1628 TGraph *gr = new TGraph(npoints);
1629 //sort xVariable increasing
1630 Int_t *sortIndex = new Int_t[npoints];
1631 TMath::Sort(npoints,x,sortIndex);
7fb602b1 1632 for (Int_t i=0;i<npoints;++i){
75d8233f 1633 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1634 }
1635
1636
1637 if ( xVariable == 2 ) delete xVar;
1638 delete x;
1639 delete y;
1640 delete sortIndex;
1641 return gr;
1642}
1643//_____________________________________________________________________
1644void AliTPCCalibCE::Analyse()
1645{
1646 //
1647 // Calculate calibration constants
1648 //
1649
1650 TVectorD paramQ(3);
1651 TVectorD paramT0(3);
1652 TVectorD paramRMS(3);
1653 TMatrixD dummy(3,3);
1654
7fb602b1 1655 for (Int_t iSec=0; iSec<72; ++iSec){
75d8233f 1656 TH2S *hT0 = GetHistoT0(iSec);
1657 if (!hT0 ) continue;
1658
1659 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1660 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1661 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1662 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1663
1664 TH2S *hQ = GetHistoQ(iSec);
1665 TH2S *hRMS = GetHistoRMS(iSec);
1666
3cd27a08 1667 Short_t *arrayhQ = hQ->GetArray();
1668 Short_t *arrayhT0 = hT0->GetArray();
1669 Short_t *arrayhRMS = hRMS->GetArray();
75d8233f 1670
1671 UInt_t nChannels = fROC->GetNChannels(iSec);
1672
1673 //debug
1674 Int_t row=0;
1675 Int_t pad=0;
1676 Int_t padc=0;
1677 //! debug
1678
7fb602b1 1679 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
75d8233f 1680
1681
1682 Float_t cogTime0 = -1000;
1683 Float_t cogQ = -1000;
1684 Float_t cogRMS = -1000;
1685 Float_t cogOut = 0;
1686
1687
1688 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1689 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1690 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1691
3cd27a08 1692 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1693 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1694 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
75d8233f 1695
1696
1697
1698 /*
7fb602b1 1699 //outlier specifications
75d8233f 1700 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1701 cogOut = 1;
1702 cogTime0 = 0;
1703 cogQ = 0;
1704 cogRMS = 0;
1705 }
1706*/
1707 rocQ->SetValue(iChannel, cogQ*cogQ);
1708 rocT0->SetValue(iChannel, cogTime0);
1709 rocRMS->SetValue(iChannel, cogRMS);
1710 rocOut->SetValue(iChannel, cogOut);
1711
1712
1713 //debug
1714 if ( fDebugLevel > 0 ){
1715 if ( !fDebugStreamer ) {
1716 //debug stream
1717 TDirectory *backup = gDirectory;
1718 fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1719 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1720 }
1721
1722 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1723 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1724 padc = pad-(fROC->GetNPads(iSec,row)/2);
1725
1726 (*fDebugStreamer) << "DataEnd" <<
1727 "Sector=" << iSec <<
1728 "Pad=" << pad <<
1729 "PadC=" << padc <<
1730 "Row=" << row <<
1731 "PadSec=" << iChannel <<
1732 "Q=" << cogQ <<
1733 "T0=" << cogTime0 <<
1734 "RMS=" << cogRMS <<
1735 "\n";
1736 }
1737 //! debug
1738
1739 }
1740
1741 }
7fb602b1 1742 if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
75d8233f 1743// delete fDebugStreamer;
bf57d87d 1744// fDebugStreamer = 0x0;
75d8233f 1745}
1746//_____________________________________________________________________
1747void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1748{
1749 //
1750 // Write class to file
1751 //
1752
1753 TString sDir(dir);
1754 TString option;
1755
1756 if ( append )
1757 option = "update";
1758 else
1759 option = "recreate";
1760
1761 TDirectory *backup = gDirectory;
1762 TFile f(filename,option.Data());
1763 f.cd();
1764 if ( !sDir.IsNull() ){
1765 f.mkdir(sDir.Data());
1766 f.cd(sDir);
1767 }
1768 this->Write();
1769 f.Close();
1770
1771 if ( backup ) backup->cd();
1772}