]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCalibSignal.cxx
Fix for negative sigmas. Memory leaks (Marco)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibSignal.cxx
CommitLineData
7e684c91 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>
7e684c91 36#include <TH2S.h>
7e684c91 37#include <TString.h>
38#include <TVectorF.h>
39#include <TMath.h>
7e684c91 40
7e684c91 41#include <TDirectory.h>
42#include <TSystem.h>
43#include <TFile.h>
44
45//AliRoot includes
46#include "AliRawReader.h"
47#include "AliRawReaderRoot.h"
bc331d5b 48#include "AliRawReaderDate.h"
7e684c91 49#include "AliTPCRawStream.h"
50#include "AliTPCCalROC.h"
51#include "AliTPCCalPad.h"
52#include "AliTPCROC.h"
53#include "AliTPCParam.h"
54#include "AliTPCCalibSignal.h"
55#include "AliTPCcalibDB.h"
bc331d5b 56#include "AliMathBase.h"
7e684c91 57#include "TTreeStream.h"
bc331d5b 58
59//date
60#include "event.h"
7e684c91 61ClassImp(AliTPCCalibSignal) /*FOLD00*/
62
63AliTPCCalibSignal::AliTPCCalibSignal() : /*FOLD00*/
64 TObject(),
65 fFirstTimeBin(60),
66 fLastTimeBin(120),
67 fFirstTimeBinT0(-15),
68 fLastTimeBinT0(15),
bc331d5b 69 fNbinsT0(200),
70 fXminT0(-2),
71 fXmaxT0(2),
72 fNbinsQ(200),
73 fXminQ(14),
74 fXmaxQ(55),
75 fNbinsRMS(100),
76 fXminRMS(0),
77 fXmaxRMS(5),
7e684c91 78 fLastSector(-1),
79 fROC(AliTPCROC::Instance()),
80 fParam(new AliTPCParam),
bc331d5b 81 fPedestalTPC(0x0),
7e684c91 82 fBpedestal(kFALSE),
83 fCalRocArrayT0(72),
84 fCalRocArrayQ(72),
85 fCalRocArrayRMS(72),
86 fCalRocArrayOutliers(72),
87 fHistoQArray(72),
88 fHistoT0Array(72),
89 fHistoRMSArray(72),
90 fPadTimesArrayEvent(72),
91 fPadQArrayEvent(72),
92 fPadRMSArrayEvent(72),
93 fPadPedestalArrayEvent(72),
94 fCurrentChannel(-1),
95 fCurrentSector(-1),
96 fCurrentRow(-1),
97 fMaxPadSignal(-1),
98 fMaxTimeBin(-1),
99 fPadSignal(1024),
100 fVTime0Offset1(72),
101 fVTime0Offset1Counter(72),
102 fEvent(-1),
103 fDebugStreamer(0x0),
104 fDebugLevel(0)
105{
106 //
107 // AliTPCSignal default constructor
108 //
109
7e684c91 110}
111//_____________________________________________________________________
112AliTPCCalibSignal::AliTPCCalibSignal(const AliTPCCalibSignal &sig) :
113 TObject(sig),
114 fFirstTimeBin(sig.fFirstTimeBin),
115 fLastTimeBin(sig.fLastTimeBin),
116 fFirstTimeBinT0(sig.fFirstTimeBinT0),
117 fLastTimeBinT0(sig.fLastTimeBinT0),
bc331d5b 118 fNbinsT0(sig.fNbinsT0),
119 fXminT0(sig.fXminT0),
120 fXmaxT0(sig.fXmaxT0),
121 fNbinsQ(sig.fNbinsQ),
122 fXminQ(sig.fXminQ),
123 fXmaxQ(sig.fXmaxQ),
124 fNbinsRMS(sig.fNbinsRMS),
125 fXminRMS(sig.fXminRMS),
126 fXmaxRMS(sig.fXmaxRMS),
7e684c91 127 fLastSector(-1),
128 fROC(AliTPCROC::Instance()),
129 fParam(new AliTPCParam),
130 fPedestalTPC(sig.fPedestalTPC),
131 fBpedestal(sig.fBpedestal),
132 fCalRocArrayT0(72),
133 fCalRocArrayQ(72),
134 fCalRocArrayRMS(72),
135 fCalRocArrayOutliers(72),
136 fHistoQArray(72),
137 fHistoT0Array(72),
138 fHistoRMSArray(72),
139 fPadTimesArrayEvent(72),
140 fPadQArrayEvent(72),
141 fPadRMSArrayEvent(72),
142 fPadPedestalArrayEvent(72),
143 fCurrentChannel(-1),
144 fCurrentSector(-1),
145 fCurrentRow(-1),
146 fMaxPadSignal(-1),
147 fMaxTimeBin(-1),
148 fPadSignal(1024),
149 fVTime0Offset1(72),
150 fVTime0Offset1Counter(72),
151 fEvent(-1),
152 fDebugStreamer(0x0),
153 fDebugLevel(sig.fDebugLevel)
154{
155 //
156 // AliTPCSignal default constructor
157 //
158
159 for (Int_t iSec = 0; iSec < 72; iSec++){
160 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
161 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
162 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
163 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
164
165 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
166 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
167 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
168
169 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
170 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
171 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
172 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
173
174 if ( hQ != 0x0 ){
175 TH2S *hNew = new TH2S(*hQ);
176 hNew->SetDirectory(0);
177 fHistoQArray.AddAt(hNew,iSec);
178 }
179 if ( hT0 != 0x0 ){
180 TH2S *hNew = new TH2S(*hT0);
181 hNew->SetDirectory(0);
182 fHistoQArray.AddAt(hNew,iSec);
183 }
184 if ( hRMS != 0x0 ){
185 TH2S *hNew = new TH2S(*hRMS);
186 hNew->SetDirectory(0);
187 fHistoQArray.AddAt(hNew,iSec);
188 }
189 }
190
7e684c91 191}
192//_____________________________________________________________________
193AliTPCCalibSignal& AliTPCCalibSignal::operator = (const AliTPCCalibSignal &source)
194{
195 //
196 // assignment operator
197 //
198 if (&source == this) return *this;
199 new (this) AliTPCCalibSignal(source);
200
201 return *this;
202}
203//_____________________________________________________________________
204AliTPCCalibSignal::~AliTPCCalibSignal()
205{
206 //
207 // destructor
208 //
209
210 fCalRocArrayT0.Delete();
211 fCalRocArrayQ.Delete();
212 fCalRocArrayRMS.Delete();
213
214 fHistoQArray.Delete();
215 fHistoT0Array.Delete();
216 fHistoRMSArray.Delete();
217
218 fPadTimesArrayEvent.Delete();
219 fPadQArrayEvent.Delete();
220 fPadRMSArrayEvent.Delete();
221 fPadPedestalArrayEvent.Delete();
222
223 if ( fDebugStreamer) delete fDebugStreamer;
224 delete fROC;
225 delete fParam;
226}
227//_____________________________________________________________________
228Int_t AliTPCCalibSignal::Update(const Int_t icsector, /*FOLD00*/
229 const Int_t icRow,
230 const Int_t icPad,
231 const Int_t icTimeBin,
232 const Float_t csignal)
233{
234 //
235 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
236 // no extra analysis necessary. Assumes knowledge of the signal shape!
237 // assumes that it is looped over consecutive time bins of one pad
238 //
239 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
240
241 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
242
243 //init first pad and sector in this event
244 if ( fCurrentChannel == -1 ) {
245 fCurrentChannel = iChannel;
246 fCurrentSector = icsector;
247 fCurrentRow = icRow;
248 }
249
250 //process last pad if we change to a new one
251 if ( iChannel != fCurrentChannel ){
252 ProcessPad();
253 fCurrentChannel = iChannel;
254 fCurrentSector = icsector;
255 fCurrentRow = icRow;
256 }
257
7e684c91 258 //fill signals for current pad
259 fPadSignal[icTimeBin]=csignal;
260 if ( csignal > fMaxPadSignal ){
261 fMaxPadSignal = csignal;
262 fMaxTimeBin = icTimeBin;
263 }
264 return 0;
265}
266//_____________________________________________________________________
267void AliTPCCalibSignal::ProcessPad() /*FOLD00*/
268{
269 //
270 // Process data of current pad
271 //
7e684c91 272
273 Float_t pedestal = 0;
274
275 if ( fBpedestal ){
276 //!!!!!!! does not work like this
277 //use pedestal database
278 AliTPCCalROC *pedestalROC = 0x0;
279
280 //only load new pedestals if the sector has changed
281 if ( fCurrentSector!=fLastSector ){
282 pedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
283 fLastSector=fCurrentSector;
284 }
285
286 pedestal = pedestalROC->GetValue(fCurrentChannel);
287
288 } else {
7e684c91 289
290 //find pedestal for pad on the fly
291 //using a few timebins before the signal
292 Int_t pminus1 = 10, pminus2=5;
293 Float_t sumN=0;
294
295 for (Int_t i=fMaxTimeBin-pminus1; i<fMaxTimeBin-pminus2+1; i++){
296 if ( i>fFirstTimeBin && i<fLastTimeBin ){
297 pedestal+=fPadSignal[i];
298 sumN+=1.;
299 }
300 }
301
302 if ( sumN>0 ) pedestal/=sumN;
7e684c91 303 }
304
305
306
307 //!!!! check borders
308 //find signal mean and sigma
309 Int_t tminus = 2, tplus=7;
310 Double_t meanT=0, sigmaT=0, Qsum=0;
311
7e684c91 312
313 for (Int_t i=fMaxTimeBin-tminus; i<fMaxTimeBin+tplus; i++){
bc331d5b 314 if ( i>=fFirstTimeBin && i<=fLastTimeBin ){
7e684c91 315 Double_t val=fPadSignal[i]-pedestal;
316 meanT+=val*(i+.5); //+.5: center of the timebin
317 sigmaT+=val*(i+.5)*(i+.5);
318 Qsum+=val;
319 }
320 }
321
322
323
324 //!!!! What to do if Qsum == 0???
325 //!!!! Should there be some threshold for max - pedestal and/or Qsum???
326 //!!!! What if Qsum < 0
327 //!!!! only fill time0 offset if Qsum > 0???
328 if ( Qsum > 0 ){
329 meanT/=Qsum;
330 sigmaT/=Qsum;
331 sigmaT = TMath::Sqrt(TMath::Abs(meanT*meanT - sigmaT));
332
333 //fill Time0 offset data for this event
334 fVTime0Offset1[fCurrentSector]+=meanT;
335 fVTime0Offset1Counter[fCurrentSector]++;
336 } else {
337 Qsum=0;
338 meanT = fLastTimeBinT0+1; //put to overflow bin
339 sigmaT = fLastTimeBinT0-fFirstTimeBinT0; //put to overflow bin
340 }
341
7e684c91 342 //Fill Event T0 counter
343 (*GetPadTimesEvent(fCurrentSector,kTRUE))[fCurrentChannel] = meanT;
344
345
346 //Normalise Q to pad area of irocs
347 Float_t norm = fParam->GetPadPitchWidth(0)*fParam->GetPadPitchLength(0,0)/(
348 fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow));
349
7e684c91 350 //Fill Q histogram
bc331d5b 351 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(Qsum*norm), fCurrentChannel );
7e684c91 352
353 //Fill RMS histogram
bc331d5b 354 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
7e684c91 355
356
357 //Fill debugging info
358 if ( fDebugLevel>0 ){
359 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=pedestal;
360 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
361 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=Qsum;
362 }
363
7e684c91 364 ResetPad();
7e684c91 365}
366//_____________________________________________________________________
367void AliTPCCalibSignal::EndEvent() /*FOLD00*/
368{
369 //
370 // Process data of current pad
371 //
bc331d5b 372 //check if last pad has allready been processed, if not do so
373 if ( fMaxTimeBin>-1 ) ProcessPad();
374
375 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
7e684c91 376 for ( Int_t iSec = 0; iSec<72; iSec++ ){
377 TVectorF *vTimes = GetPadTimesEvent(iSec);
378 if ( !vTimes ) continue;
379
380 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); iChannel++ ){
381 Float_t Time0 = fVTime0Offset1[iSec]/fVTime0Offset1Counter[iSec];
382 Float_t Time = (*vTimes)[iChannel];
383
bc331d5b 384 GetHistoT0(iSec,kTRUE)->Fill( Time-Time0,iChannel );
7e684c91 385
386
387 //Debug start
388 if ( fDebugLevel>0 ){
bc331d5b 389 if ( !fDebugStreamer ) {
390 //debug stream
391 TDirectory *backup = gDirectory;
392 fDebugStreamer = new TTreeSRedirector("deb2.root");
393 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
394 }
395
7e684c91 396 Int_t row=0;
397 Int_t pad=0;
398 Int_t padc=0;
399
400 Float_t Q = (*GetPadQEvent(iSec))[iChannel];
401 Float_t RMS = (*GetPadRMSEvent(iSec))[iChannel];
402
403 UInt_t channel=iChannel;
404 Int_t sector=iSec;
405
406 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
407 pad = channel-fROC->GetRowIndexes(sector)[row];
408 padc = pad-(fROC->GetNPads(sector,row)/2);
409
410 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
411 Form("hSignalD%d.%d.%d",sector,row,pad),
412 fLastTimeBin-fFirstTimeBin,
413 fFirstTimeBin,fLastTimeBin);
414 h1->SetDirectory(0);
415
416 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
417 h1->Fill(i,fPadSignal(i));
418
419 (*fDebugStreamer) << "DataPad" <<
420 "Event=" << fEvent <<
421 "Sector="<< sector <<
422 "Row=" << row<<
423 "Pad=" << pad <<
424 "PadC=" << padc <<
425 "PadSec="<< channel <<
426 "Time0=" << Time0 <<
427 "Time=" << Time <<
428 "RMS=" << RMS <<
429 "Sum=" << Q <<
430 "hist.=" << h1 <<
431 "\n";
432
433 delete h1;
434 }
435 //Debug end
436
437 }
438 }
439
440}
441//_____________________________________________________________________
bc331d5b 442Bool_t AliTPCCalibSignal::ProcessEvent(AliTPCRawStream *rawStream) /*FOLD00*/
7e684c91 443{
444 //
bc331d5b 445 // Event Processing loop - AliTPCRawStream
7e684c91 446 //
7e684c91 447
bc331d5b 448 rawStream->SetOldRCUFormat(1);
7e684c91 449
450 ResetEvent();
451
452 Bool_t withInput = kFALSE;
453
bc331d5b 454 while (rawStream->Next()) {
7e684c91 455
bc331d5b 456 Int_t isector = rawStream->GetSector(); // current sector
457 Int_t iRow = rawStream->GetRow(); // current row
458 Int_t iPad = rawStream->GetPad(); // current pad
459 Int_t iTimeBin = rawStream->GetTime(); // current time bin
460 Float_t signal = rawStream->GetSignal(); // current ADC signal
7e684c91 461
462 Update(isector,iRow,iPad,iTimeBin,signal);
463 withInput = kTRUE;
464 }
465
466 if (withInput){
7e684c91 467 EndEvent();
468 }
469
7e684c91 470 return withInput;
471}
472//_____________________________________________________________________
bc331d5b 473Bool_t AliTPCCalibSignal::ProcessEvent(AliRawReader *rawReader)
474{
475 //
476 // Event processing loop - AliRawReader
477 //
478
479
480 AliTPCRawStream rawStream(rawReader);
481
482 rawReader->Select("TPC");
483
484 return ProcessEvent(&rawStream);
485}
486//_____________________________________________________________________
487Bool_t AliTPCCalibSignal::ProcessEvent(eventHeaderStruct *event)
488{
489 //
490 // Event processing loop - date event
491 //
492 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
493 Bool_t result=ProcessEvent(rawReader);
494 delete rawReader;
495 return result;
496
497}
498//_____________________________________________________________________
7e684c91 499TH2S* AliTPCCalibSignal::GetHisto(Int_t sector, TObjArray *arr, /*FOLD00*/
500 Int_t nbinsY, Float_t ymin, Float_t ymax,
501 Char_t *type, Bool_t force)
502{
503 //
504 // return pointer to Q histogram
505 // if force is true create a new histogram if it doesn't exist allready
506 //
507 if ( !force || arr->UncheckedAt(sector) )
508 return (TH2S*)arr->UncheckedAt(sector);
509
510 // if we are forced and histogram doesn't yes exist create it
511 Char_t name[255], title[255];
512
513 sprintf(name,"hCalib%s%.2d",type,sector);
514 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
515
516 // new histogram with Q calib information. One value for each pad!
517 TH2S* hist = new TH2S(name,title,
bc331d5b 518 nbinsY, ymin, ymax,
519 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
7e684c91 520 hist->SetDirectory(0);
521 arr->AddAt(hist,sector);
522 return hist;
523}
524//_____________________________________________________________________
525TH2S* AliTPCCalibSignal::GetHistoT0(Int_t sector, Bool_t force) /*FOLD00*/
526{
527 //
528 // return pointer to T0 histogram
529 // if force is true create a new histogram if it doesn't exist allready
530 //
531 TObjArray *arr = &fHistoT0Array;
bc331d5b 532 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
7e684c91 533}
534//_____________________________________________________________________
535TH2S* AliTPCCalibSignal::GetHistoQ(Int_t sector, Bool_t force) /*FOLD00*/
536{
537 //
538 // return pointer to Q histogram
539 // if force is true create a new histogram if it doesn't exist allready
540 //
541 TObjArray *arr = &fHistoQArray;
bc331d5b 542 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
7e684c91 543}
544//_____________________________________________________________________
545TH2S* AliTPCCalibSignal::GetHistoRMS(Int_t sector, Bool_t force) /*FOLD00*/
546{
547 //
548 // return pointer to Q histogram
549 // if force is true create a new histogram if it doesn't exist allready
550 //
551 TObjArray *arr = &fHistoRMSArray;
bc331d5b 552 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
7e684c91 553}
554//_____________________________________________________________________
555TVectorF* AliTPCCalibSignal::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force) /*FOLD00*/
556{
557 //
558 // return pointer to Pad Info from 'arr' for the current event and sector
559 // if force is true create it if it doesn't exist allready
560 //
561 if ( !force || arr->UncheckedAt(sector) )
562 return (TVectorF*)arr->UncheckedAt(sector);
563
564 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
565 arr->AddAt(vect,sector);
566 return vect;
567}
568//_____________________________________________________________________
569TVectorF* AliTPCCalibSignal::GetPadTimesEvent(Int_t sector, Bool_t force) /*FOLD00*/
570{
571 //
572 // return pointer to Pad Times Array for the current event and sector
573 // if force is true create it if it doesn't exist allready
574 //
575 TObjArray *arr = &fPadTimesArrayEvent;
576 return GetPadInfoEvent(sector,arr,force);
577}
578//_____________________________________________________________________
579TVectorF* AliTPCCalibSignal::GetPadQEvent(Int_t sector, Bool_t force) /*FOLD00*/
580{
581 //
582 // return pointer to Pad Q Array for the current event and sector
583 // if force is true create it if it doesn't exist allready
584 // for debugging purposes only
585 //
586
587 TObjArray *arr = &fPadQArrayEvent;
588 return GetPadInfoEvent(sector,arr,force);
589}
590//_____________________________________________________________________
591TVectorF* AliTPCCalibSignal::GetPadRMSEvent(Int_t sector, Bool_t force) /*FOLD00*/
592{
593 //
594 // return pointer to Pad RMS Array for the current event and sector
595 // if force is true create it if it doesn't exist allready
596 // for debugging purposes only
597 //
598 TObjArray *arr = &fPadRMSArrayEvent;
599 return GetPadInfoEvent(sector,arr,force);
600}
601//_____________________________________________________________________
602TVectorF* AliTPCCalibSignal::GetPadPedestalEvent(Int_t sector, Bool_t force) /*FOLD00*/
603{
604 //
605 // return pointer to Pad RMS Array for the current event and sector
606 // if force is true create it if it doesn't exist allready
607 // for debugging purposes only
608 //
609 TObjArray *arr = &fPadPedestalArrayEvent;
610 return GetPadInfoEvent(sector,arr,force);
611}
612//_____________________________________________________________________
613AliTPCCalROC* AliTPCCalibSignal::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) /*FOLD00*/
614{
615 //
616 // return pointer to ROC Calibration
617 // if force is true create a new histogram if it doesn't exist allready
618 //
619 if ( !force || arr->UncheckedAt(sector) )
620 return (AliTPCCalROC*)arr->UncheckedAt(sector);
621
622 // if we are forced and histogram doesn't yes exist create it
623
624 // new AliTPCCalROC for T0 information. One value for each pad!
625 AliTPCCalROC *croc = new AliTPCCalROC(sector);
626 //init values
627 for ( UInt_t iChannel = 0; iChannel<croc->GetNchannels(); iChannel++){
628 croc->SetValue(iChannel, 0);
629 }
630 arr->AddAt(croc,sector);
631 return croc;
632}
633//_____________________________________________________________________
634AliTPCCalROC* AliTPCCalibSignal::GetCalRocT0(Int_t sector, Bool_t force) /*FOLD00*/
635{
636 //
637 // return pointer to Carge ROC Calibration
638 // if force is true create a new histogram if it doesn't exist allready
639 //
640 TObjArray *arr = &fCalRocArrayT0;
641 return GetCalRoc(sector, arr, force);
642}
643//_____________________________________________________________________
644AliTPCCalROC* AliTPCCalibSignal::GetCalRocQ(Int_t sector, Bool_t force) /*FOLD00*/
645{
646 //
647 // return pointer to T0 ROC Calibration
648 // if force is true create a new histogram if it doesn't exist allready
649 //
650 TObjArray *arr = &fCalRocArrayQ;
651 return GetCalRoc(sector, arr, force);
652}
653//_____________________________________________________________________
654AliTPCCalROC* AliTPCCalibSignal::GetCalRocRMS(Int_t sector, Bool_t force) /*FOLD00*/
655{
656 //
657 // return pointer to signal width ROC Calibration
658 // if force is true create a new histogram if it doesn't exist allready
659 //
660 TObjArray *arr = &fCalRocArrayRMS;
661 return GetCalRoc(sector, arr, force);
662}
663//_____________________________________________________________________
664AliTPCCalROC* AliTPCCalibSignal::GetCalRocOutliers(Int_t sector, Bool_t force)
665{
666 //
667 // return pointer to Outliers
668 // if force is true create a new histogram if it doesn't exist allready
669 //
670 TObjArray *arr = &fCalRocArrayOutliers;
671 return GetCalRoc(sector, arr, force);
672}
673//_____________________________________________________________________
674void AliTPCCalibSignal::ResetEvent() /*FOLD00*/
675{
676 //
677 // Reset global counters -- Should be called before each event is processed
678 //
679 fLastSector=-1;
680 fCurrentSector=-1;
681 fCurrentRow=-1;
682 fCurrentChannel=-1;
683
684 ResetPad();
685
686 fPadTimesArrayEvent.Delete();
687 fPadQArrayEvent.Delete();
688 fPadRMSArrayEvent.Delete();
689 fPadPedestalArrayEvent.Delete();
690
691 for ( Int_t i=0; i<72; i++ ){
692 fVTime0Offset1[i]=0;
693 fVTime0Offset1Counter[i]=0;
694 }
695}
696//_____________________________________________________________________
697void AliTPCCalibSignal::ResetPad() /*FOLD00*/
698{
699 //
700 // Reset pad infos -- Should be called after a pad has been processed
701 //
702 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; i++)
703 fPadSignal[i] = 0;
704 fMaxTimeBin = -1;
705 fMaxPadSignal = -1;
706}
707//_____________________________________________________________________
708void AliTPCCalibSignal::Analyse() /*FOLD00*/
709{
710 //
711 // Calculate calibration constants
712 //
713
bc331d5b 714 TVectorD paramQ(3);
715 TVectorD paramT0(3);
716 TVectorD paramRMS(3);
717 TMatrixD dummy(3,3);
7e684c91 718
719 for (Int_t iSec=0; iSec<72; iSec++){
720 TH2S *hT0 = GetHistoT0(iSec);
721 if (!hT0 ) continue;
722
723 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
724 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
725 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
726 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
727
728 TH2S *hQ = GetHistoQ(iSec);
729 TH2S *hRMS = GetHistoRMS(iSec);
730
bc331d5b 731 Short_t *array_hQ = hQ->GetArray();
732 Short_t *array_hT0 = hT0->GetArray();
733 Short_t *array_hRMS = hRMS->GetArray();
734
735 UInt_t nChannels = fROC->GetNChannels(iSec);
736
7e684c91 737 //debug
738 Int_t row=0;
739 Int_t pad=0;
740 Int_t padc=0;
741 //! debug
742
bc331d5b 743 for (UInt_t iChannel=0; iChannel<nChannels; iChannel++){
7e684c91 744
745
746 Float_t cogTime0 = -1000;
747 Float_t cogQ = -1000;
748 Float_t cogRMS = -1000;
749 Float_t cogOut = 0;
750
7e684c91 751
bc331d5b 752 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
753 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
754 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
755
7e684c91 756/*
bc331d5b 757 AliMathBase::FitGaus(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
758 AliMathBase::FitGaus(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
759 AliMathBase::FitGaus(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
760 cogQ = paramQ[1];
761 cogTime0 = paramT0[1];
762 cogRMS = paramRMS[1];
763*/
764 cogQ = AliMathBase::GetCOG(array_hQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
765 cogTime0 = AliMathBase::GetCOG(array_hT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
766 cogRMS = AliMathBase::GetCOG(array_hRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
767
768
769
770 /*
7e684c91 771 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
772 cogOut = 1;
773 cogTime0 = 0;
774 cogQ = 0;
775 cogRMS = 0;
776 }
777*/
778 rocQ->SetValue(iChannel, cogQ*cogQ);
779 rocT0->SetValue(iChannel, cogTime0);
780 rocRMS->SetValue(iChannel, cogRMS);
781 rocOut->SetValue(iChannel, cogOut);
782
783
784 //debug
785 if ( fDebugLevel > 0 ){
bc331d5b 786 if ( !fDebugStreamer ) {
787 //debug stream
788 TDirectory *backup = gDirectory;
789 fDebugStreamer = new TTreeSRedirector("deb2.root");
790 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
791 }
792
7e684c91 793 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
794 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
795 padc = pad-(fROC->GetNPads(iSec,row)/2);
796
797 (*fDebugStreamer) << "DataEnd" <<
798 "Sector=" << iSec <<
799 "Pad=" << pad <<
800 "PadC=" << padc <<
801 "Row=" << row <<
802 "PadSec=" << iChannel <<
803 "Q=" << cogQ <<
804 "T0=" << cogTime0 <<
805 "RMS=" << cogRMS <<
806 "\n";
807 }
808 //! debug
809
810 }
811
812 }
813 delete fDebugStreamer;
814 fDebugStreamer = 0x0;
815}
816//_____________________________________________________________________
817void AliTPCCalibSignal::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
818{
819 //
820 // Write class to file
821 //
822
7e684c91 823 TString sDir(dir);
824 TString option;
825
826 if ( append )
827 option = "update";
828 else
829 option = "recreate";
830
bc331d5b 831 TDirectory *backup = gDirectory;
7e684c91 832 TFile f(filename,option.Data());
bc331d5b 833 f.cd();
7e684c91 834 if ( !sDir.IsNull() ){
835 f.mkdir(sDir.Data());
836 f.cd(sDir);
837 }
bc331d5b 838 this->Write();
7e684c91 839 f.Close();
840
841 if ( backup ) backup->cd();
7e684c91 842}