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