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