]> 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
Data base generation class for configuration entries (Haavard)
[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
16
17
18
19
20//-------------------------------------------------------
21// Implementation of the TPC pulser calibration
22//
23// Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch
24//
25//
26//-------------------------------------------------------
27
28
29/* $Id$ */
30
31
32
33//Root includes
34#include <TObjArray.h>
35#include <TH1F.h>
36#include <TH2S.h>
37#include <TString.h>
38#include <TVectorF.h>
39#include <TVectorD.h>
40#include <TMatrixD.h>
41#include <TMath.h>
42#include <TGraph.h>
43
44#include <TDirectory.h>
45#include <TSystem.h>
46#include <TFile.h>
47
48//AliRoot includes
49#include "AliRawReader.h"
50#include "AliRawReaderRoot.h"
51#include "AliRawReaderDate.h"
52#include "AliRawEventHeaderBase.h"
53#include "AliTPCRawStream.h"
54#include "AliTPCcalibDB.h"
55#include "AliTPCCalROC.h"
56#include "AliTPCCalPad.h"
57#include "AliTPCROC.h"
58#include "AliTPCParam.h"
59#include "AliTPCCalibCE.h"
75d8233f 60#include "AliMathBase.h"
61#include "TTreeStream.h"
62
63//date
64#include "event.h"
65ClassImp(AliTPCCalibCE) /*FOLD00*/
66
67AliTPCCalibCE::AliTPCCalibCE() : /*FOLD00*/
68 TObject(),
69 fFirstTimeBin(650),
70 fLastTimeBin(1000),
71 fNbinsT0(200),
72 fXminT0(-5),
73 fXmaxT0(5),
74 fNbinsQ(200),
bf57d87d 75 fXminQ(1),
76 fXmaxQ(40),
75d8233f 77 fNbinsRMS(100),
bf57d87d 78 fXminRMS(0.1),
79 fXmaxRMS(5.1),
75d8233f 80 fLastSector(-1),
81 fOldRCUformat(kTRUE),
82 fROC(AliTPCROC::Instance()),
83 fParam(new AliTPCParam),
84 fPedestalTPC(0x0),
85 fPadNoiseTPC(0x0),
86 fPedestalROC(0x0),
87 fPadNoiseROC(0x0),
75d8233f 88 fCalRocArrayT0(72),
89 fCalRocArrayQ(72),
90 fCalRocArrayRMS(72),
91 fCalRocArrayOutliers(72),
92 fHistoQArray(72),
93 fHistoT0Array(72),
94 fHistoRMSArray(72),
95 fHistoTmean(72),
96 fParamArrayEvent(1000),
97 fParamArrayEventPol1(72),
98 fParamArrayEventPol2(72),
99 fTMeanArrayEvent(1000),
bf57d87d 100 fQMeanArrayEvent(1000),
75d8233f 101 fVEventTime(1000),
102 fVEventNumber(1000),
103 fNevents(0),
104 fTimeStamp(0),
105 fRunNumber(-1),
106 fOldRunNumber(-1),
107 fPadTimesArrayEvent(72),
108 fPadQArrayEvent(72),
109 fPadRMSArrayEvent(72),
110 fPadPedestalArrayEvent(72),
111 fCurrentChannel(-1),
112 fCurrentSector(-1),
113 fCurrentRow(-1),
114 fMaxPadSignal(-1),
115 fMaxTimeBin(-1),
116 fPadSignal(1024),
117 fPadPedestal(0),
118 fPadNoise(0),
119 fVTime0Offset(72),
120 fVTime0OffsetCounter(72),
bf57d87d 121 fVMeanQ(72),
122 fVMeanQCounter(72),
75d8233f 123// fHTime0(0x0),
124 fEvent(-1),
125 fDebugStreamer(0x0),
126 fDebugLevel(0)
127{
128 //
129 // AliTPCSignal default constructor
130 //
131// fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
132}
133//_____________________________________________________________________
134AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
135 TObject(sig),
136 fFirstTimeBin(sig.fFirstTimeBin),
137 fLastTimeBin(sig.fLastTimeBin),
138 fNbinsT0(sig.fNbinsT0),
139 fXminT0(sig.fXminT0),
140 fXmaxT0(sig.fXmaxT0),
141 fNbinsQ(sig.fNbinsQ),
142 fXminQ(sig.fXminQ),
143 fXmaxQ(sig.fXmaxQ),
144 fNbinsRMS(sig.fNbinsRMS),
145 fXminRMS(sig.fXminRMS),
146 fXmaxRMS(sig.fXmaxRMS),
147 fLastSector(-1),
148 fOldRCUformat(kTRUE),
149 fROC(AliTPCROC::Instance()),
150 fParam(new AliTPCParam),
151 fPedestalTPC(0x0),
152 fPadNoiseTPC(0x0),
153 fPedestalROC(0x0),
154 fPadNoiseROC(0x0),
75d8233f 155 fCalRocArrayT0(72),
156 fCalRocArrayQ(72),
157 fCalRocArrayRMS(72),
158 fCalRocArrayOutliers(72),
159 fHistoQArray(72),
160 fHistoT0Array(72),
161 fHistoRMSArray(72),
162 fHistoTmean(72),
163 fParamArrayEvent(1000),
164 fParamArrayEventPol1(72),
165 fParamArrayEventPol2(72),
166 fTMeanArrayEvent(1000),
bf57d87d 167 fQMeanArrayEvent(1000),
75d8233f 168 fVEventTime(1000),
169 fVEventNumber(1000),
170 fNevents(sig.fNevents),
171 fTimeStamp(0),
172 fRunNumber(-1),
173 fOldRunNumber(-1),
174 fPadTimesArrayEvent(72),
175 fPadQArrayEvent(72),
176 fPadRMSArrayEvent(72),
177 fPadPedestalArrayEvent(72),
178 fCurrentChannel(-1),
179 fCurrentSector(-1),
180 fCurrentRow(-1),
181 fMaxPadSignal(-1),
182 fMaxTimeBin(-1),
183 fPadSignal(1024),
184 fPadPedestal(0),
bf57d87d 185 fPadNoise(0),
75d8233f 186 fVTime0Offset(72),
187 fVTime0OffsetCounter(72),
bf57d87d 188 fVMeanQ(72),
189 fVMeanQCounter(72),
75d8233f 190// fHTime0(0x0),
191 fEvent(-1),
192 fDebugStreamer(0x0),
193 fDebugLevel(sig.fDebugLevel)
194{
195 //
196 // AliTPCSignal default constructor
197 //
198
199 for (Int_t iSec = 0; iSec < 72; iSec++){
200 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
201 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
202 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
203 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
204
205 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
206 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
207 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
208
209 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
210 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
211 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
212 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
213
214 if ( hQ != 0x0 ){
215 TH2S *hNew = new TH2S(*hQ);
216 hNew->SetDirectory(0);
217 fHistoQArray.AddAt(hNew,iSec);
218 }
219 if ( hT0 != 0x0 ){
220 TH2S *hNew = new TH2S(*hT0);
221 hNew->SetDirectory(0);
222 fHistoQArray.AddAt(hNew,iSec);
223 }
224 if ( hRMS != 0x0 ){
225 TH2S *hNew = new TH2S(*hRMS);
226 hNew->SetDirectory(0);
227 fHistoQArray.AddAt(hNew,iSec);
228 }
229 }
230 for (Int_t iEvent=0; iEvent<sig.fParamArrayEvent.GetSize(); iEvent++)
231 fParamArrayEvent.AddAtAndExpand(sig.fParamArrayEvent.At(iEvent),iEvent);
232 Int_t nrows = sig.fVEventTime.GetNrows();
233 fVEventTime.ResizeTo(nrows);
234 for (Int_t iEvent=0; iEvent<nrows; iEvent++)
235 fVEventTime[iEvent] = sig.fVEventTime[iEvent];
236
237// fHTime0 = new TH1F("hTime0Event","hTime0Event",(fLastTimeBin-fFirstTimeBin)*10,fFirstTimeBin,fLastTimeBin);
238}
239//_____________________________________________________________________
240AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
241{
242 //
243 // assignment operator
244 //
245 if (&source == this) return *this;
246 new (this) AliTPCCalibCE(source);
247
248 return *this;
249}
250//_____________________________________________________________________
251AliTPCCalibCE::~AliTPCCalibCE()
252{
253 //
254 // destructor
255 //
256
257 fCalRocArrayT0.Delete();
258 fCalRocArrayQ.Delete();
259 fCalRocArrayRMS.Delete();
260
261 fHistoQArray.Delete();
262 fHistoT0Array.Delete();
263 fHistoRMSArray.Delete();
264
265 fPadTimesArrayEvent.Delete();
266 fPadQArrayEvent.Delete();
267 fPadRMSArrayEvent.Delete();
268 fPadPedestalArrayEvent.Delete();
269
270 if ( fDebugStreamer) delete fDebugStreamer;
271// if ( fHTime0 ) delete fHTime0;
272 delete fROC;
273 delete fParam;
274}
275//_____________________________________________________________________
276Int_t AliTPCCalibCE::Update(const Int_t icsector, /*FOLD00*/
277 const Int_t icRow,
278 const Int_t icPad,
279 const Int_t icTimeBin,
280 const Float_t csignal)
281{
282 //
283 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
284 // no extra analysis necessary. Assumes knowledge of the signal shape!
285 // assumes that it is looped over consecutive time bins of one pad
286 //
287 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
288
289 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
290
291 //init first pad and sector in this event
292 if ( fCurrentChannel == -1 ) {
293 fCurrentChannel = iChannel;
294 fCurrentSector = icsector;
295 fCurrentRow = icRow;
296 }
297
298 //process last pad if we change to a new one
299 if ( iChannel != fCurrentChannel ){
300 ProcessPad();
301 fCurrentChannel = iChannel;
302 fCurrentSector = icsector;
303 fCurrentRow = icRow;
304 }
305
306 //fill signals for current pad
bf57d87d 307 fPadSignal.GetMatrixArray()[icTimeBin]=csignal;
75d8233f 308 if ( csignal > fMaxPadSignal ){
309 fMaxPadSignal = csignal;
310 fMaxTimeBin = icTimeBin;
311 }
312 return 0;
313}
314//_____________________________________________________________________
315void AliTPCCalibCE::FindPedestal(Float_t part)
316{
317 //
318 // find pedestal and noise for the current pad. Use either database or
319 // truncated mean with part*100%
320 //
bf57d87d 321 Bool_t noPedestal = kTRUE;;
322 if (fPedestalTPC&&fPadNoiseTPC){
75d8233f 323 //use pedestal database
75d8233f 324 //only load new pedestals if the sector has changed
325 if ( fCurrentSector!=fLastSector ){
326 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
327 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
328 fLastSector=fCurrentSector;
329 }
330
bf57d87d 331 if ( fPedestalROC&&fPadNoiseROC ){
75d8233f 332 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
333 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
334 noPedestal = kFALSE;
335 }
336
337 }
338
339 //if we are not running with pedestal database, or for the current sector there is no information
340 //available, calculate the pedestal and noise on the fly
341 if ( noPedestal ) {
342 const Int_t kPedMax = 100; //maximum pedestal value
343 Float_t max = 0;
344 Float_t maxPos = 0;
345 Int_t median = -1;
346 Int_t count0 = 0;
347 Int_t count1 = 0;
348 //
bf57d87d 349 Float_t padSignal=0;
350 //
75d8233f 351 UShort_t histo[kPedMax];
352 memset(histo,0,kPedMax*sizeof(UShort_t));
bf57d87d 353
75d8233f 354 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; i++){
bf57d87d 355 padSignal = fPadSignal.GetMatrixArray()[i];
356 if (padSignal<=0) continue;
357 if (padSignal>max && i>10) {
358 max = padSignal;
75d8233f 359 maxPos = i;
360 }
bf57d87d 361 if (padSignal>kPedMax-1) continue;
362 histo[int(padSignal+0.5)]++;
75d8233f 363 count0++;
364 }
365 //
366 for (Int_t i=1; i<kPedMax; i++){
367 if (count1<count0*0.5) median=i;
368 count1+=histo[i];
369 }
370 // truncated mean
371 //
372 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
373 //
374 for (Int_t idelta=1; idelta<10; idelta++){
375 if (median-idelta<=0) continue;
376 if (median+idelta>kPedMax) continue;
377 if (count<part*count1){
378 count+=histo[median-idelta];
379 mean +=histo[median-idelta]*(median-idelta);
380 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
381 count+=histo[median+idelta];
382 mean +=histo[median+idelta]*(median+idelta);
383 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
384 }
385 }
386 fPadPedestal = 0;
387 fPadNoise = 0;
388 if ( count > 0 ) {
389 mean/=count;
390 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
391 fPadPedestal = mean;
392 fPadNoise = rms;
393 }
394 }
395}
396//_____________________________________________________________________
397void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
398{
399 //
400 // Find position, signal width and height of the CE signal (last signal)
401 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
402 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
403 //
404
405 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
406 Int_t cemaxpos = 0;
407 Float_t ceSumThreshold = 8.*fPadNoise; // threshold for the signal sum
408 const Int_t kCemin = 4; // range for the analysis of the ce signal +- channels from the peak
409 const Int_t kCemax = 7;
410
bf57d87d 411 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
412
75d8233f 413 // find maximum closest to the sector mean from the last event
414 for ( Int_t imax=0; imax<maxima.GetNrows(); imax++){
bf57d87d 415 Float_t tmean = (*((TVectorF*)(fTMeanArrayEvent[fTMeanArrayEvent.GetLast()])))[fCurrentSector];
75d8233f 416 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
417 minDist = tmean-maxima[imax];
418 cemaxpos = (Int_t)maxima[imax];
419 }
420 }
421
422 if (cemaxpos!=0){
bf57d87d 423 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
75d8233f 424 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; i++){
bf57d87d 425 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
426 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
427 ceTime+=signal*(i+0.5);
428 ceRMS +=signal*(i+0.5)*(i+0.5);
429 ceQsum+=signal;
75d8233f 430 }
431 }
432 }
433 if (ceQmax&&ceQsum>ceSumThreshold) {
434 ceTime/=ceQsum;
435 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
bf57d87d 436 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
437 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
438
439 //Normalise Q to pad area of irocs
440 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
441
442 ceQsum/=norm;
443 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
444 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
75d8233f 445 } else {
446 ceQmax=0;
447 ceTime=0;
448 ceRMS =0;
449 ceQsum=0;
450 }
451 param[0] = ceQmax;
452 param[1] = ceTime;
453 param[2] = ceRMS;
454 qSum = ceQsum;
455}
456//_____________________________________________________________________
457Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus)
458{
459 //
460 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
461 // and 'tplus' timebins after 'pos'
462 //
463 for (Int_t iTime = pos; iTime>pos-tminus; iTime--)
464 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
465 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; iTime++,iTime2++){
466 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
467 iTime2++;
468 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
469 }
470 return kTRUE;
471}
472//_____________________________________________________________________
473void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
474{
475 //
476 // Find local maxima on the pad signal and Histogram them
477 //
478 Float_t ceThreshold = 5.*fPadNoise; // threshold for the signal
479 Int_t count = 0;
480 Int_t tminus = 2;
481 Int_t tplus = 3;
482 for (Int_t i=fLastTimeBin-tplus-1; i>=fFirstTimeBin+tminus; i--){
483 if ( (fPadSignal[i]-fPadPedestal)>ceThreshold && IsPeak(i,tminus,tplus) ){
bf57d87d 484 maxima.GetMatrixArray()[count++]=i;
75d8233f 485 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
486 }
487 }
488}
489//_____________________________________________________________________
490void AliTPCCalibCE::ProcessPad() /*FOLD00*/
491{
492 //
493 // Process data of current pad
494 //
495 FindPedestal();
496
497 TVectorF maxima(10);
498 FindLocalMaxima(maxima);
499 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
500
501
502
503 TVectorD param(3);
504 Float_t Qsum;
505 FindCESignal(param, Qsum, maxima);
506
507 Double_t meanT = param[1];
508 Double_t sigmaT = param[2];
509
510 //Fill Event T0 counter
bf57d87d 511 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
75d8233f 512
513 //Fill Q histogram
bf57d87d 514 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum), fCurrentChannel );
75d8233f 515
516 //Fill RMS histogram
517 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
518
519
520 //Fill debugging info
521 if ( fDebugLevel>0 ){
bf57d87d 522 (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
523 (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
524 (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=Qsum;
75d8233f 525 }
526
527 ResetPad();
528}
529//_____________________________________________________________________
530void AliTPCCalibCE::EndEvent() /*FOLD00*/
531{
532 //
533 // Process data of current pad
534 // The Functions 'SetTimeStamp' and 'SetRunNumber' should be called
535 // before the EndEvent function to set the event timestamp and number!!!
536 // This is automatically done if the ProcessEvent(AliRawReader *rawReader)
537 // function was called
538 //
539
540 //check if last pad has allready been processed, if not do so
541 if ( fMaxTimeBin>-1 ) ProcessPad();
542
543 TVectorD param(3);
544 TMatrixD dummy(3,3);
545 TVectorF vMeanTime(72);
bf57d87d 546 TVectorF vMeanQ(72);
75d8233f 547 AliTPCCalROC calIroc(0);
548 AliTPCCalROC calOroc(36);
549
550 //find mean time0 offset for side A and C
551 Double_t time0Side[2]; //time0 for side A:0 and C:0
552 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:0
553 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
554 for ( Int_t iSec = 0; iSec<72; iSec++ ){
bf57d87d 555 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
556 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
75d8233f 557 }
bf57d87d 558 if ( time0SideCount[0] >0 )
559 time0Side[0]/=time0SideCount[0];
560 if ( time0SideCount[1] >0 )
561 time0Side[1]/=time0SideCount[1];
75d8233f 562 // end find time0 offset
bf57d87d 563
75d8233f 564 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
565 for ( Int_t iSec = 0; iSec<72; iSec++ ){
75d8233f 566
bf57d87d 567 //find median and then calculate the mean around it
568 TH1S *hMeanT = GetHistoTmean(iSec);
569 if ( !hMeanT ) continue;
75d8233f 570 Double_t entries = hMeanT->GetEntries();
571 Double_t sum = 0;
572 Short_t *arr = hMeanT->GetArray()+1;
573 Int_t ibin=0;
574 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ibin++){
575 sum+=arr[ibin];
576 if ( sum>=(entries/2) ) break;
577 }
578 Int_t delta = 4;
579 Int_t firstBin = fFirstTimeBin+ibin-delta;
580 Int_t lastBin = fFirstTimeBin+ibin+delta;
581 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
582 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
bf57d87d 583 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
584 vMeanTime.GetMatrixArray()[iSec]=median;
585 // end find median
586
587 TVectorF *vTimes = GetPadTimesEvent(iSec);
588 if ( !vTimes ) continue;
589 AliTPCCalROC calIrocOutliers(0);
590 AliTPCCalROC calOrocOutliers(36);
591
592 // calculate mean Q of the sector
593 Float_t meanQ = 0;
594 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
595 vMeanQ.GetMatrixArray()[iSec]=meanQ;
75d8233f 596
597 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++ ){
bf57d87d 598 Float_t Time = (*vTimes).GetMatrixArray()[iChannel];
75d8233f 599
600 //set values for temporary roc calibration class
bf57d87d 601 if ( iSec < 36 ) {
602 calIroc.SetValue(iChannel, Time);
603 if ( Time == 0 ) calIrocOutliers.SetValue(iChannel,1);
604
605 } else {
606 calOroc.SetValue(iChannel, Time);
607 if ( Time == 0 ) calOrocOutliers.SetValue(iChannel,1);
608 }
75d8233f 609
610 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
611 GetHistoT0(iSec,kTRUE)->Fill( Time-time0Side[(iSec/18)%2],iChannel );
612
613
614
615 //------------------------------- Debug start ------------------------------
616 if ( fDebugLevel>0 ){
617 if ( !fDebugStreamer ) {
618 //debug stream
619 TDirectory *backup = gDirectory;
620 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
621 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
622 }
623
624 Int_t row=0;
625 Int_t pad=0;
626 Int_t padc=0;
627
628 Float_t Q = (*GetPadQEvent(iSec))[iChannel];
629 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
630
631 UInt_t channel=iChannel;
632 Int_t sector=iSec;
633
634 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
635 pad = channel-fROC->GetRowIndexes(sector)[row];
636 padc = pad-(fROC->GetNPads(sector,row)/2);
637
638// TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
639// Form("hSignalD%d.%d.%d",sector,row,pad),
640// fLastTimeBin-fFirstTimeBin,
641// fFirstTimeBin,fLastTimeBin);
642// h1->SetDirectory(0);
643//
644// for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
645// h1->Fill(i,fPadSignal(i));
646
bf57d87d 647 Double_t T0Sec = 0;
648 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
649 T0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
650 Double_t T0Side = time0Side[(iSec/18)%2];
75d8233f 651 (*fDebugStreamer) << "DataPad" <<
652 "Event=" << fNevents <<
653 "EventNumber=" << fRunNumber <<
bf57d87d 654 "TimeStamp=" << fTimeStamp <<
75d8233f 655 "Sector="<< sector <<
656 "Row=" << row<<
657 "Pad=" << pad <<
658 "PadC=" << padc <<
659 "PadSec="<< channel <<
660 "Time0Sec=" << T0Sec <<
bf57d87d 661 "Time0Side=" << T0Side <<
75d8233f 662 "Time=" << Time <<
663 "RMS=" << RMS <<
664 "Sum=" << Q <<
bf57d87d 665 "MeanQ=" << meanQ <<
666 // "hist.=" << h1 <<
75d8233f 667 "\n";
668
bf57d87d 669 // delete h1;
670
75d8233f 671 }
672 //----------------------------- Debug end ------------------------------
673 }// end channel loop
674 hMeanT->Reset();
675
676 TVectorD paramPol1(3);
677 TVectorD paramPol2(6);
678 TMatrixD matPol1(3,3);
679 TMatrixD matPol2(6,6);
bf57d87d 680 Float_t chi2Pol1=0;
681 Float_t chi2Pol2=0;
75d8233f 682
683 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
684 if ( iSec < 36 ){
bf57d87d 685 calIroc.GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
686 calIroc.GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
75d8233f 687 } else {
bf57d87d 688 calOroc.GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
689 calOroc.GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
75d8233f 690 }
691
692 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
693 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
694 }
bf57d87d 695// printf("events: %d -- size: %d\n",fNevents,GetParamArrayPol1(iSec)->GetSize());
75d8233f 696
697 //------------------------------- Debug start ------------------------------
698 if ( fDebugLevel>0 ){
bf57d87d 699 if ( !fDebugStreamer ) {
700 //debug stream
701 TDirectory *backup = gDirectory;
702 fDebugStreamer = new TTreeSRedirector("debugCalibCE.root");
703 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
704 }
75d8233f 705 (*fDebugStreamer) << "DataRoc" <<
706 "Event=" << fEvent <<
707 "EventNumber=" << fRunNumber <<
708 "TimeStamp=" << fTimeStamp <<
709 "Sector="<< iSec <<
710 "hMeanT.=" << hMeanT <<
711 "median=" << median <<
712 "paramPol1.=" << &paramPol1 <<
713 "paramPol2.=" << &paramPol2 <<
714 "matPol1.=" << &matPol1 <<
715 "matPol2.=" << &matPol2 <<
716 "chi2Pol1=" << chi2Pol1 <<
717 "chi2Pol2=" << chi2Pol2 <<
718 "\n";
719 }
720 //------------------------------- Debug end ------------------------------
721 }// end sector loop
722
723/* AliMathBase::FitGaus(fHTime0->GetArray()+1,
724 fHTime0->GetNbinsX(),
725 fHTime0->GetXaxis()->GetXmin(),
726 fHTime0->GetXaxis()->GetXmax(),
727 &param, &dummy);*/
728// fHTime0->Reset();
729
730 // fParamArrayEvent.AddAtAndExpand(new TVectorD(param),fNevents);
731 fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
bf57d87d 732 fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
75d8233f 733 if ( fVEventTime.GetNrows() < fNevents ) {
734 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+1000));
735 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+1000));
736 }
737 fVEventTime[fNevents] = fTimeStamp;
738 fVEventNumber[fNevents] = fRunNumber;
739
740 fNevents++;
741 fOldRunNumber = fRunNumber;
742
743}
744//_____________________________________________________________________
745Bool_t AliTPCCalibCE::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
746{
747 //
748 // Event Processing loop - AliTPCRawStream
749 // The Function 'SetTimeStamp' should be called for each event to set the event time stamp!!!
750 //
751
752 rawStream->SetOldRCUFormat(fOldRCUformat);
753
754 ResetEvent();
755
756 Bool_t withInput = kFALSE;
757
758 while (rawStream->Next()) {
759
760 Int_t isector = rawStream->GetSector(); // current sector
761 Int_t iRow = rawStream->GetRow(); // current row
762 Int_t iPad = rawStream->GetPad(); // current pad
763 Int_t iTimeBin = rawStream->GetTime(); // current time bin
764 Float_t signal = rawStream->GetSignal(); // current ADC signal
765
766 Update(isector,iRow,iPad,iTimeBin,signal);
767 withInput = kTRUE;
768 }
769
770 if (withInput){
771 EndEvent();
772 }
773
774 return withInput;
775}
776//_____________________________________________________________________
777Bool_t AliTPCCalibCE::ProcessEvent(AliRawReader *rawReader)
778{
779 //
780 // Event processing loop - AliRawReader
781 //
782
783
784 AliTPCRawStream rawStream(rawReader);
785 AliRawEventHeaderBase* eventHeader = (AliRawEventHeaderBase*)rawReader->GetEventHeader();
786 if (eventHeader){
787 fTimeStamp = eventHeader->Get("Timestamp");
788 fRunNumber = eventHeader->Get("RunNb");
789 }
790
791
792 rawReader->Select("TPC");
793
794 return ProcessEvent(&rawStream);
795}
796//_____________________________________________________________________
797Bool_t AliTPCCalibCE::ProcessEvent(eventHeaderStruct *event)
798{
799 //
800 // Event processing loop - date event
801 //
802 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
803 Bool_t result=ProcessEvent(rawReader);
804 delete rawReader;
805 return result;
806
807}
808//_____________________________________________________________________
809TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
810 Int_t nbinsY, Float_t ymin, Float_t ymax,
811 Char_t *type, Bool_t force)
812{
813 //
814 // return pointer to TH2S histogram of 'type'
815 // if force is true create a new histogram if it doesn't exist allready
816 //
817 if ( !force || arr->UncheckedAt(sector) )
818 return (TH2S*)arr->UncheckedAt(sector);
819
820 // if we are forced and histogram doesn't yes exist create it
821 Char_t name[255], title[255];
822
823 sprintf(name,"hCalib%s%.2d",type,sector);
824 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
825
826 // new histogram with Q calib information. One value for each pad!
827 TH2S* hist = new TH2S(name,title,
828 nbinsY, ymin, ymax,
829 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
830 hist->SetDirectory(0);
831 arr->AddAt(hist,sector);
832 return hist;
833}
834//_____________________________________________________________________
835TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
836{
837 //
838 // return pointer to T0 histogram
839 // if force is true create a new histogram if it doesn't exist allready
840 //
841 TObjArray *arr = &fHistoT0Array;
842 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
843}
844//_____________________________________________________________________
845TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
846{
847 //
848 // return pointer to Q histogram
849 // if force is true create a new histogram if it doesn't exist allready
850 //
851 TObjArray *arr = &fHistoQArray;
852 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
853}
854//_____________________________________________________________________
855TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
856{
857 //
858 // return pointer to Q histogram
859 // if force is true create a new histogram if it doesn't exist allready
860 //
861 TObjArray *arr = &fHistoRMSArray;
862 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
863}
864//_____________________________________________________________________
865TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
866 Char_t *type, Bool_t force)
867{
868 //
869 // return pointer to TH1S histogram
870 // if force is true create a new histogram if it doesn't exist allready
871 //
872 if ( !force || arr->UncheckedAt(sector) )
873 return (TH1S*)arr->UncheckedAt(sector);
874
875 // if we are forced and histogram doesn't yes exist create it
876 Char_t name[255], title[255];
877
878 sprintf(name,"hCalib%s%.2d",type,sector);
879 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
880
881 // new histogram with Q calib information. One value for each pad!
882 TH1S* hist = new TH1S(name,title,
883 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
884 hist->SetDirectory(0);
885 arr->AddAt(hist,sector);
886 return hist;
887}
888//_____________________________________________________________________
889TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
890{
891 //
892 // return pointer to Q histogram
893 // if force is true create a new histogram if it doesn't exist allready
894 //
895 TObjArray *arr = &fHistoTmean;
896 return GetHisto(sector, arr, "LastTmean", force);
897}
898//_____________________________________________________________________
899TVectorF* AliTPCCalibCE::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
900{
901 //
902 // return pointer to Pad Info from 'arr' for the current event and sector
903 // if force is true create it if it doesn't exist allready
904 //
905 if ( !force || arr->UncheckedAt(sector) )
906 return (TVectorF*)arr->UncheckedAt(sector);
907
908 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
909 arr->AddAt(vect,sector);
910 return vect;
911}
912//_____________________________________________________________________
913TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
914{
915 //
916 // return pointer to Pad Times Array for the current event and sector
917 // if force is true create it if it doesn't exist allready
918 //
919 TObjArray *arr = &fPadTimesArrayEvent;
920 return GetPadInfoEvent(sector,arr,force);
921}
922//_____________________________________________________________________
923TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
924{
925 //
926 // return pointer to Pad Q Array for the current event and sector
927 // if force is true create it if it doesn't exist allready
928 // for debugging purposes only
929 //
930
931 TObjArray *arr = &fPadQArrayEvent;
932 return GetPadInfoEvent(sector,arr,force);
933}
934//_____________________________________________________________________
935TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
936{
937 //
938 // return pointer to Pad RMS Array for the current event and sector
939 // if force is true create it if it doesn't exist allready
940 // for debugging purposes only
941 //
942 TObjArray *arr = &fPadRMSArrayEvent;
943 return GetPadInfoEvent(sector,arr,force);
944}
945//_____________________________________________________________________
946TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
947{
948 //
949 // return pointer to Pad RMS Array for the current event and sector
950 // if force is true create it if it doesn't exist allready
951 // for debugging purposes only
952 //
953 TObjArray *arr = &fPadPedestalArrayEvent;
954 return GetPadInfoEvent(sector,arr,force);
955}
956//_____________________________________________________________________
957AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
958{
959 //
960 // return pointer to ROC Calibration
961 // if force is true create a new histogram if it doesn't exist allready
962 //
963 if ( !force || arr->UncheckedAt(sector) )
964 return (AliTPCCalROC*)arr->UncheckedAt(sector);
965
966 // if we are forced and histogram doesn't yes exist create it
967
968 // new AliTPCCalROC for T0 information. One value for each pad!
969 AliTPCCalROC *croc = new AliTPCCalROC(sector);
970 //init values
971 for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
972 croc->SetValue(iChannel, 0);
973 }
974 arr->AddAt(croc,sector);
975 return croc;
976}
977//_____________________________________________________________________
978AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
979{
980 //
981 // return pointer to Carge ROC Calibration
982 // if force is true create a new histogram if it doesn't exist allready
983 //
984 TObjArray *arr = &fCalRocArrayT0;
985 return GetCalRoc(sector, arr, force);
986}
987//_____________________________________________________________________
988AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
989{
990 //
991 // return pointer to T0 ROC Calibration
992 // if force is true create a new histogram if it doesn't exist allready
993 //
994 TObjArray *arr = &fCalRocArrayQ;
995 return GetCalRoc(sector, arr, force);
996}
997//_____________________________________________________________________
998AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
999{
1000 //
1001 // return pointer to signal width ROC Calibration
1002 // if force is true create a new histogram if it doesn't exist allready
1003 //
1004 TObjArray *arr = &fCalRocArrayRMS;
1005 return GetCalRoc(sector, arr, force);
1006}
1007//_____________________________________________________________________
1008AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1009{
1010 //
1011 // return pointer to Outliers
1012 // if force is true create a new histogram if it doesn't exist allready
1013 //
1014 TObjArray *arr = &fCalRocArrayOutliers;
1015 return GetCalRoc(sector, arr, force);
1016}
1017//_____________________________________________________________________
1018TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force)
1019{
1020 //
1021 // return pointer to TObjArray of fit parameters
1022 // if force is true create a new histogram if it doesn't exist allready
1023 //
1024 if ( !force || arr->UncheckedAt(sector) )
1025 return (TObjArray*)arr->UncheckedAt(sector);
1026
1027 // if we are forced and array doesn't yes exist create it
1028
1029 // new TObjArray for parameters
1030 TObjArray *newArr = new TObjArray;
1031 arr->AddAt(newArr,sector);
1032 return newArr;
1033}
1034//_____________________________________________________________________
1035TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1036{
1037 //
1038 // return pointer to TObjArray of fit parameters from plane fit
1039 // if force is true create a new histogram if it doesn't exist allready
1040 //
1041 TObjArray *arr = &fParamArrayEventPol1;
1042 return GetParamArray(sector, arr, force);
1043}
1044//_____________________________________________________________________
1045TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1046{
1047 //
1048 // return pointer to TObjArray of fit parameters from parabola fit
1049 // if force is true create a new histogram if it doesn't exist allready
1050 //
1051 TObjArray *arr = &fParamArrayEventPol2;
1052 return GetParamArray(sector, arr, force);
1053}
1054//_____________________________________________________________________
1055void AliTPCCalibCE::ResetEvent() /*FOLD00*/
1056{
1057 //
1058 // Reset global counters -- Should be called before each event is processed
1059 //
1060 fLastSector=-1;
1061 fCurrentSector=-1;
1062 fCurrentRow=-1;
1063 fCurrentChannel=-1;
1064
1065 ResetPad();
1066
1067 fPadTimesArrayEvent.Delete();
1068 fPadQArrayEvent.Delete();
1069 fPadRMSArrayEvent.Delete();
1070 fPadPedestalArrayEvent.Delete();
1071
1072 for ( Int_t i=0; i<72; i++ ){
bf57d87d 1073 fVTime0Offset.GetMatrixArray()[i]=0;
1074 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1075 fVMeanQ.GetMatrixArray()[i]=0;
1076 fVMeanQCounter.GetMatrixArray()[i]=0;
75d8233f 1077 }
1078}
1079//_____________________________________________________________________
1080void AliTPCCalibCE::ResetPad() /*FOLD00*/
1081{
1082 //
1083 // Reset pad infos -- Should be called after a pad has been processed
1084 //
1085 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
bf57d87d 1086 fPadSignal.GetMatrixArray()[i] = 0;
75d8233f 1087 fMaxTimeBin = -1;
1088 fMaxPadSignal = -1;
1089 fPadPedestal = -1;
1090 fPadNoise = -1;
1091}
1092//_____________________________________________________________________
1093TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter) /*FOLD00*/
1094{
1095 //
1096 // Make graph from fit parameters of pol1 or pol2 fit
1097 // xVariable: 0-run time, 1-run number, 2-internal event counter
bf57d87d 1098 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 2-mean Q
75d8233f 1099 // fitParameter: fit parameter ( 0-2 for pol1, 0-5 for pol2, 0 for mean time )
1100 //
1101
1102 Double_t *x = new Double_t[fNevents];
1103 Double_t *y = new Double_t[fNevents];
1104
1105 TVectorD *xVar = 0x0;
1106 TObjArray *aType = 0x0;
1107 Int_t npoints=0;
1108
1109 // sanity checks
1110 if ( (sector<0) || (sector>71) ) return 0x0;
1111 if ( (xVariable<0) || (xVariable>2) ) return 0x0;
bf57d87d 1112 if ( (fitType<0) || (fitType>3) ) return 0x0;
75d8233f 1113 if ( fitType==0 ){
1114 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1115 aType = &fParamArrayEventPol1;
1116 if ( aType->At(sector)==0x0 ) return 0x0;
1117 }
1118 else if ( fitType==1 ){
1119 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1120 aType = &fParamArrayEventPol2;
1121 if ( aType->At(sector)==0x0 ) return 0x0;
1122 }
1123
1124
1125 if ( xVariable == 0 ) xVar = &fVEventTime;
1126 if ( xVariable == 1 ) xVar = &fVEventNumber;
1127 if ( xVariable == 2 ) {
1128 xVar = new TVectorD(fNevents);
1129 for ( Int_t i=0;i<fNevents; i++) (*xVar)[i]=i;
1130 }
1131
1132 for (Int_t ievent =0; ievent<fNevents; ievent++){
1133 if ( fitType<2 ){
1134 TObjArray *events = (TObjArray*)(aType->At(sector));
1135 if ( events->GetSize()<=ievent ) break;
1136 TVectorD *v = (TVectorD*)(events->At(ievent));
75d8233f 1137 if ( v!=0x0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1138 } else if (fitType == 2) {
1139 Double_t yValue=(*((TVectorF*)(fTMeanArrayEvent[ievent])))[sector];
bf57d87d 1140 if ( yValue>0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=yValue;npoints++;}
1141 }else if (fitType == 3) {
1142 Double_t yValue=(*((TVectorF*)(fQMeanArrayEvent[ievent])))[sector];
75d8233f 1143 if ( yValue>0 ) { x[npoints]=(*xVar)[ievent]; y[npoints]=yValue;npoints++;}
1144 }
1145 }
1146
1147 TGraph *gr = new TGraph(npoints);
1148 //sort xVariable increasing
1149 Int_t *sortIndex = new Int_t[npoints];
1150 TMath::Sort(npoints,x,sortIndex);
1151 for (Int_t i=0;i<npoints;i++){
1152 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1153 }
1154
1155
1156 if ( xVariable == 2 ) delete xVar;
1157 delete x;
1158 delete y;
1159 delete sortIndex;
1160 return gr;
1161}
1162//_____________________________________________________________________
1163void AliTPCCalibCE::Analyse()
1164{
1165 //
1166 // Calculate calibration constants
1167 //
1168
1169 TVectorD paramQ(3);
1170 TVectorD paramT0(3);
1171 TVectorD paramRMS(3);
1172 TMatrixD dummy(3,3);
1173
1174 for (Int_t iSec=0; iSec<72; iSec++){
1175 TH2S *hT0 = GetHistoT0(iSec);
1176 if (!hT0 ) continue;
1177
1178 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1179 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1180 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1181 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1182
1183 TH2S *hQ = GetHistoQ(iSec);
1184 TH2S *hRMS = GetHistoRMS(iSec);
1185
1186 Short_t *array_hQ = hQ->GetArray();
1187 Short_t *array_hT0 = hT0->GetArray();
1188 Short_t *array_hRMS = hRMS->GetArray();
1189
1190 UInt_t nChannels = fROC->GetNChannels(iSec);
1191
1192 //debug
1193 Int_t row=0;
1194 Int_t pad=0;
1195 Int_t padc=0;
1196 //! debug
1197
1198 for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
1199
1200
1201 Float_t cogTime0 = -1000;
1202 Float_t cogQ = -1000;
1203 Float_t cogRMS = -1000;
1204 Float_t cogOut = 0;
1205
1206
1207 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1208 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1209 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1210
1211/*
1212 AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1213 AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1214 AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
1215 cogQ = paramQ[1];
1216 cogTime0 = paramT0[1];
1217 cogRMS = paramRMS[1];
1218*/
1219 cogQ = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1220 cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1221 cogRMS = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
1222
1223
1224
1225 /*
1226 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1227 cogOut = 1;
1228 cogTime0 = 0;
1229 cogQ = 0;
1230 cogRMS = 0;
1231 }
1232*/
1233 rocQ->SetValue(iChannel, cogQ*cogQ);
1234 rocT0->SetValue(iChannel, cogTime0);
1235 rocRMS->SetValue(iChannel, cogRMS);
1236 rocOut->SetValue(iChannel, cogOut);
1237
1238
1239 //debug
1240 if ( fDebugLevel > 0 ){
1241 if ( !fDebugStreamer ) {
1242 //debug stream
1243 TDirectory *backup = gDirectory;
1244 fDebugStreamer = new TTreeSRedirector("debugCalibCEAnalysis.root");
1245 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1246 }
1247
1248 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1249 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1250 padc = pad-(fROC->GetNPads(iSec,row)/2);
1251
1252 (*fDebugStreamer) << "DataEnd" <<
1253 "Sector=" << iSec <<
1254 "Pad=" << pad <<
1255 "PadC=" << padc <<
1256 "Row=" << row <<
1257 "PadSec=" << iChannel <<
1258 "Q=" << cogQ <<
1259 "T0=" << cogTime0 <<
1260 "RMS=" << cogRMS <<
1261 "\n";
1262 }
1263 //! debug
1264
1265 }
1266
1267 }
bf57d87d 1268 fDebugStreamer->GetFile()->Write();
75d8233f 1269// delete fDebugStreamer;
bf57d87d 1270// fDebugStreamer = 0x0;
75d8233f 1271}
1272//_____________________________________________________________________
1273void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1274{
1275 //
1276 // Write class to file
1277 //
1278
1279 TString sDir(dir);
1280 TString option;
1281
1282 if ( append )
1283 option = "update";
1284 else
1285 option = "recreate";
1286
1287 TDirectory *backup = gDirectory;
1288 TFile f(filename,option.Data());
1289 f.cd();
1290 if ( !sDir.IsNull() ){
1291 f.mkdir(sDir.Data());
1292 f.cd(sDir);
1293 }
1294 this->Write();
1295 f.Close();
1296
1297 if ( backup ) backup->cd();
1298}