]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCalibPulser.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Merge Jens and Christian changes
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibPulser.cxx
CommitLineData
467a4b25 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
3cd27a08 18/////////////////////////////////////////////////////////////////////////////////////////
19// //
20// Implementation of the TPC pulser calibration //
21// //
22// Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
23// //
24/////////////////////////////////////////////////////////////////////////////////////////
467a4b25 25/***************************************************************************
26 * Class Description *
27 ***************************************************************************
28
29 The AliTPCCalibPulser class is used to get calibration data concerning the FEE using
30 runs performed with the calibration pulser.
31
32 The information retrieved is
33 - Time0 differences
34 - Signal width differences
35 - Amplification variations
36
37 the seen differences arise from the manufacturing tolerances of the PASAs and are very small within
38 one chip and somewhat large between different chips.
39
40 Histograms:
41 For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum]) is created when
42 it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are stored in the
43 TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.
44
45
46 Working principle:
47 ------------------
48 Raw calibration pulser data is processed by calling one of the ProcessEvent(...) functions
49 (see below). These in the end call the Update(...) function.
50
51 - the Update(...) function:
52 In this function the array fPadSignal is filled with the adc signals between the specified range
53 fFirstTimeBin and fLastTimeBin for the current pad.
54 before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
55 stored in fPadSignal.
56
57 - the ProcessPad() function:
58 Find Pedestal and Noise information
59 - use database information which has to be set by calling
60 SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)
61 - if no information from the pedestal data base
62 is available the informaion is calculated on the fly ( see FindPedestal() function )
63
64 Find the Pulser signal information
65 - calculate mean = T0, RMS = signal width and Q sum in a range of -2+7 timebins around Q max
66 the Q sum is scaled by pad area
67 (see FindPulserSignal(...) function)
68
69 Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)
70 Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE)),
71
72 At the end of each event the EndEvent() function is called
73
74 - the EndEvent() function:
75 calculate the mean T0 for each ROC and fill the Time0 histogram with Time0-<Time0 for ROC>
76 This is done to overcome syncronisation problems between the trigger and the fec clock.
77
78 After accumulating the desired statistics the Analyse() function has to be called.
79 - the Analyse() function
80 Whithin this function the mean values of T0, RMS, Q are calculated for each pad, using
81 the AliMathBase::GetCOG(...) function, and the calibration
82 storage classes (AliTPCCalROC) are filled for each ROC.
83 The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and
84 fCalRocArrayQ;
85
86
87
88 User interface for filling data:
89 --------------------------------
90
91 To Fill information one of the following functions can be used:
92
93 Bool_t ProcessEvent(eventHeaderStruct *event);
94 - process Date event
95 - use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)
96
97 Bool_t ProcessEvent(AliRawReader *rawReader);
98 - process AliRawReader event
99 - use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)
100
101 Bool_t ProcessEvent(AliTPCRawStream *rawStream);
102 - process event from AliTPCRawStream
103 - call Update function for signal filling
104
105 Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
106 iPad, const Int_t iTimeBin, const Float_t signal);
107 - directly fill signal information (sector, row, pad, time bin, pad)
108 to the reference histograms
109
110 It is also possible to merge two independently taken calibrations using the function
111
112 void Merge(AliTPCCalibPulser *sig)
113 - copy histograms in 'sig' if the do not exist in this instance
114 - Add histograms in 'sig' to the histograms in this instance if the allready exist
115 - After merging call Analyse again!
116
117
118
119 -- example: filling data using root raw data:
120 void fillSignal(Char_t *filename)
121 {
122 rawReader = new AliRawReaderRoot(fileName);
123 if ( !rawReader ) return;
124 AliTPCCalibPulser *calib = new AliTPCCalibPulser;
125 while (rawReader->NextEvent()){
126 calib->ProcessEvent(rawReader);
127 }
128 calib->Analyse();
129 calib->DumpToFile("SignalData.root");
130 delete rawReader;
131 delete calib;
132 }
133
134
135 What kind of information is stored and how to retrieve them:
136 ------------------------------------------------------------
137
138 - Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):
139
140 TH2F *GetHistoT0(Int_t sector);
141 TH2F *GetHistoRMS(Int_t sector);
142 TH2F *GetHistoQ(Int_t sector);
143
144 - Accessing the calibration storage objects:
145
146 AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values
147 AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values
148 AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values
149
150 example for visualisation:
151 if the file "SignalData.root" was created using the above example one could do the following:
152
153 TFile fileSignal("SignalData.root")
154 AliTPCCalibPulser *sig = (AliTPCCalibPulser*)fileSignal->Get("AliTPCCalibPulser");
155 sig->GetCalRocT0(0)->Draw("colz");
156 sig->GetCalRocRMS(0)->Draw("colz");
157
158 or use the AliTPCCalPad functionality:
159 AliTPCCalPad padT0(ped->GetCalPadT0());
160 AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
161 padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
162 padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
163*/
164
165
3cd27a08 166//Root includes
167#include <TObjArray.h>
168#include <TH1F.h>
169#include <TH2S.h>
170#include <TString.h>
171#include <TVectorF.h>
172#include <TMath.h>
173
174#include <TDirectory.h>
175#include <TSystem.h>
176#include <TROOT.h>
177#include <TFile.h>
178
179//AliRoot includes
180#include "AliRawReader.h"
181#include "AliRawReaderRoot.h"
182#include "AliRawReaderDate.h"
183#include "AliTPCRawStream.h"
08205ed7 184#include "AliTPCRawStreamFast.h"
3cd27a08 185#include "AliTPCCalROC.h"
186#include "AliTPCCalPad.h"
187#include "AliTPCROC.h"
188#include "AliTPCParam.h"
189#include "AliTPCCalibPulser.h"
190#include "AliTPCcalibDB.h"
191#include "AliMathBase.h"
192#include "TTreeStream.h"
467a4b25 193
3cd27a08 194//date
195#include "event.h"
196
197
198
199
200ClassImp(AliTPCCalibPulser)
201
202AliTPCCalibPulser::AliTPCCalibPulser() :
467a4b25 203 TObject(),
204 fFirstTimeBin(60),
205 fLastTimeBin(120),
206 fNbinsT0(200),
207 fXminT0(-2),
208 fXmaxT0(2),
209 fNbinsQ(200),
210 fXminQ(1),
211 fXmaxQ(40),
212 fNbinsRMS(100),
213 fXminRMS(0.1),
214 fXmaxRMS(5.1),
215 fLastSector(-1),
216 fOldRCUformat(kTRUE),
217 fROC(AliTPCROC::Instance()),
218 fParam(new AliTPCParam),
219 fPedestalTPC(0x0),
220 fPadNoiseTPC(0x0),
221 fPedestalROC(0x0),
222 fPadNoiseROC(0x0),
223 fCalRocArrayT0(72),
224 fCalRocArrayQ(72),
225 fCalRocArrayRMS(72),
226 fCalRocArrayOutliers(72),
227 fHistoQArray(72),
228 fHistoT0Array(72),
229 fHistoRMSArray(72),
230 fPadTimesArrayEvent(72),
231 fPadQArrayEvent(72),
232 fPadRMSArrayEvent(72),
233 fPadPedestalArrayEvent(72),
234 fCurrentChannel(-1),
235 fCurrentSector(-1),
236 fCurrentRow(-1),
237 fMaxPadSignal(-1),
238 fMaxTimeBin(-1),
239 fPadSignal(1024),
240 fPadPedestal(0),
241 fPadNoise(0),
242 fVTime0Offset(72),
243 fVTime0OffsetCounter(72),
3cd27a08 244// fEvent(-1),
467a4b25 245 fDebugStreamer(0x0),
246 fDebugLevel(0)
247{
248 //
249 // AliTPCSignal default constructor
250 //
251
252}
253//_____________________________________________________________________
254AliTPCCalibPulser::AliTPCCalibPulser(const AliTPCCalibPulser &sig) :
255 TObject(sig),
256 fFirstTimeBin(sig.fFirstTimeBin),
257 fLastTimeBin(sig.fLastTimeBin),
258 fNbinsT0(sig.fNbinsT0),
259 fXminT0(sig.fXminT0),
260 fXmaxT0(sig.fXmaxT0),
261 fNbinsQ(sig.fNbinsQ),
262 fXminQ(sig.fXminQ),
263 fXmaxQ(sig.fXmaxQ),
264 fNbinsRMS(sig.fNbinsRMS),
265 fXminRMS(sig.fXminRMS),
266 fXmaxRMS(sig.fXmaxRMS),
267 fLastSector(-1),
268 fOldRCUformat(kTRUE),
269 fROC(AliTPCROC::Instance()),
270 fParam(new AliTPCParam),
271 fPedestalTPC(0x0),
272 fPadNoiseTPC(0x0),
273 fPedestalROC(0x0),
274 fPadNoiseROC(0x0),
275 fCalRocArrayT0(72),
276 fCalRocArrayQ(72),
277 fCalRocArrayRMS(72),
278 fCalRocArrayOutliers(72),
279 fHistoQArray(72),
280 fHistoT0Array(72),
281 fHistoRMSArray(72),
282 fPadTimesArrayEvent(72),
283 fPadQArrayEvent(72),
284 fPadRMSArrayEvent(72),
285 fPadPedestalArrayEvent(72),
286 fCurrentChannel(-1),
287 fCurrentSector(-1),
288 fCurrentRow(-1),
289 fMaxPadSignal(-1),
290 fMaxTimeBin(-1),
291 fPadSignal(1024),
292 fPadPedestal(0),
293 fPadNoise(0),
294 fVTime0Offset(72),
295 fVTime0OffsetCounter(72),
3cd27a08 296// fEvent(-1),
467a4b25 297 fDebugStreamer(0x0),
298 fDebugLevel(sig.fDebugLevel)
299{
300 //
301 // AliTPCSignal default constructor
302 //
303
304 for (Int_t iSec = 0; iSec < 72; ++iSec){
305 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
306 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
307 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
308 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
309
310 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
311 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
312 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
313
314 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
315 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
316 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
317 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
318
319 if ( hQ != 0x0 ){
320 TH2S *hNew = new TH2S(*hQ);
321 hNew->SetDirectory(0);
322 fHistoQArray.AddAt(hNew,iSec);
323 }
324 if ( hT0 != 0x0 ){
325 TH2S *hNew = new TH2S(*hT0);
326 hNew->SetDirectory(0);
327 fHistoQArray.AddAt(hNew,iSec);
328 }
329 if ( hRMS != 0x0 ){
330 TH2S *hNew = new TH2S(*hRMS);
331 hNew->SetDirectory(0);
332 fHistoQArray.AddAt(hNew,iSec);
333 }
334 }
335
336}
337//_____________________________________________________________________
338AliTPCCalibPulser& AliTPCCalibPulser::operator = (const AliTPCCalibPulser &source)
339{
340 //
341 // assignment operator
342 //
343 if (&source == this) return *this;
344 new (this) AliTPCCalibPulser(source);
345
346 return *this;
347}
348//_____________________________________________________________________
349AliTPCCalibPulser::~AliTPCCalibPulser()
350{
351 //
352 // destructor
353 //
354
355 fCalRocArrayT0.Delete();
356 fCalRocArrayQ.Delete();
357 fCalRocArrayRMS.Delete();
358 fCalRocArrayOutliers.Delete();
359
360 fHistoQArray.Delete();
361 fHistoT0Array.Delete();
362 fHistoRMSArray.Delete();
363
364 fPadTimesArrayEvent.Delete();
365 fPadQArrayEvent.Delete();
366 fPadRMSArrayEvent.Delete();
367 fPadPedestalArrayEvent.Delete();
368
369 if ( fDebugStreamer) delete fDebugStreamer;
3cd27a08 370 delete fROC;
467a4b25 371 delete fParam;
372}
373//_____________________________________________________________________
3cd27a08 374Int_t AliTPCCalibPulser::Update(const Int_t icsector,
467a4b25 375 const Int_t icRow,
376 const Int_t icPad,
377 const Int_t icTimeBin,
378 const Float_t csignal)
379{
380 //
3cd27a08 381 // Signal filling methode on the fly pedestal and time offset correction if necessary.
467a4b25 382 // no extra analysis necessary. Assumes knowledge of the signal shape!
383 // assumes that it is looped over consecutive time bins of one pad
384 //
385 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
386
387 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
388
389 //init first pad and sector in this event
390 if ( fCurrentChannel == -1 ) {
391 fCurrentChannel = iChannel;
392 fCurrentSector = icsector;
393 fCurrentRow = icRow;
394 }
395
396 //process last pad if we change to a new one
397 if ( iChannel != fCurrentChannel ){
398 ProcessPad();
399 fCurrentChannel = iChannel;
400 fCurrentSector = icsector;
401 fCurrentRow = icRow;
402 }
403
404 //fill signals for current pad
405 fPadSignal[icTimeBin]=csignal;
406 if ( csignal > fMaxPadSignal ){
407 fMaxPadSignal = csignal;
408 fMaxTimeBin = icTimeBin;
409 }
410 return 0;
411}
412//_____________________________________________________________________
413void AliTPCCalibPulser::FindPedestal(Float_t part)
414{
415 //
416 // find pedestal and noise for the current pad. Use either database or
417 // truncated mean with part*100%
418 //
419 Bool_t noPedestal = kTRUE;;
420 if (fPedestalTPC&&fPadNoiseTPC){
421 //use pedestal database
422 //only load new pedestals if the sector has changed
423 if ( fCurrentSector!=fLastSector ){
424 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
425 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
426 fLastSector=fCurrentSector;
427 }
428
429 if ( fPedestalROC&&fPadNoiseROC ){
430 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel);
431 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
432 noPedestal = kFALSE;
433 }
434
435 }
436
437 //if we are not running with pedestal database, or for the current sector there is no information
438 //available, calculate the pedestal and noise on the fly
439 if ( noPedestal ) {
440 const Int_t kPedMax = 100; //maximum pedestal value
441 Float_t max = 0;
442 Float_t maxPos = 0;
443 Int_t median = -1;
444 Int_t count0 = 0;
445 Int_t count1 = 0;
446 //
447 Float_t padSignal=0;
448 //
449 UShort_t histo[kPedMax];
450 memset(histo,0,kPedMax*sizeof(UShort_t));
451
452 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
453 padSignal = fPadSignal.GetMatrixArray()[i];
454 if (padSignal<=0) continue;
455 if (padSignal>max && i>10) {
456 max = padSignal;
457 maxPos = i;
458 }
459 if (padSignal>kPedMax-1) continue;
460 histo[Int_t(padSignal+0.5)]++;
461 count0++;
462 }
463 //
464 for (Int_t i=1; i<kPedMax; ++i){
465 if (count1<count0*0.5) median=i;
466 count1+=histo[i];
467 }
468 // truncated mean
469 //
470 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
471 //
472 for (Int_t idelta=1; idelta<10; ++idelta){
473 if (median-idelta<=0) continue;
474 if (median+idelta>kPedMax) continue;
475 if (count<part*count1){
476 count+=histo[median-idelta];
477 mean +=histo[median-idelta]*(median-idelta);
478 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
479 count+=histo[median+idelta];
480 mean +=histo[median+idelta]*(median+idelta);
481 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
482 }
483 }
484 fPadPedestal = 0;
485 fPadNoise = 0;
486 if ( count > 0 ) {
487 mean/=count;
488 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
489 fPadPedestal = mean;
490 fPadNoise = rms;
491 }
492 }
493}
494//_____________________________________________________________________
495void AliTPCCalibPulser::FindPulserSignal(TVectorD &param, Float_t &qSum)
496{
497 //
498 // Find position, signal width and height of the CE signal (last signal)
499 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
500 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
501 //
502
503 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
504 Int_t cemaxpos = fMaxTimeBin;
505 Float_t ceSumThreshold = 8.*fPadNoise; // threshold for the signal sum
506 const Int_t kCemin = 2; // range for the analysis of the ce signal +- channels from the peak
507 const Int_t kCemax = 7;
508
509 if (cemaxpos!=0){
510 ceQmax = fPadSignal.GetMatrixArray()[cemaxpos]-fPadPedestal;
511 for (Int_t i=cemaxpos-kCemin; i<cemaxpos+kCemax; ++i){
512 Float_t signal = fPadSignal.GetMatrixArray()[i]-fPadPedestal;
513 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) && (signal>0) ){
514 ceTime+=signal*(i+0.5);
515 ceRMS +=signal*(i+0.5)*(i+0.5);
516 ceQsum+=signal;
517 }
518 }
519 }
520 if (ceQmax&&ceQsum>ceSumThreshold) {
521 ceTime/=ceQsum;
522 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
523 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
524 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
525
526 //Normalise Q to the pad area
527 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
528
529 ceQsum/=norm;
530 } else {
531 ceQmax=0;
532 ceTime=0;
533 ceRMS =0;
534 ceQsum=0;
535 }
536 param[0] = ceQmax;
537 param[1] = ceTime;
538 param[2] = ceRMS;
539 qSum = ceQsum;
540}
541//_____________________________________________________________________
3cd27a08 542void AliTPCCalibPulser::ProcessPad()
467a4b25 543{
544 //
545 // Process data of current pad
546 //
547
548 FindPedestal();
549 TVectorD param(3);
3cd27a08 550 Float_t qSum;
551 FindPulserSignal(param, qSum);
467a4b25 552
553 Double_t meanT = param[1];
554 Double_t sigmaT = param[2];
555
556 //Fill Event T0 counter
557 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
558
559 //Fill Q histogram
3cd27a08 560 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
467a4b25 561
562 //Fill RMS histogram
563 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
564
565
566 //Fill debugging info
567 if ( fDebugLevel>0 ){
568 (*GetPadPedestalEvent(fCurrentSector,kTRUE))[fCurrentChannel]=fPadPedestal;
569 (*GetPadRMSEvent(fCurrentSector,kTRUE))[fCurrentChannel]=sigmaT;
3cd27a08 570 (*GetPadQEvent(fCurrentSector,kTRUE))[fCurrentChannel]=qSum;
467a4b25 571 }
467a4b25 572 ResetPad();
573}
574//_____________________________________________________________________
3cd27a08 575void AliTPCCalibPulser::EndEvent()
467a4b25 576{
577 //
578 // Process data of current event
579 //
580
581 //check if last pad has allready been processed, if not do so
582 if ( fMaxTimeBin>-1 ) ProcessPad();
583
584 //loop over all ROCs, fill Time0 histogram corrected for the mean Time0 of each ROC
585 for ( Int_t iSec = 0; iSec<72; ++iSec ){
586 TVectorF *vTimes = GetPadTimesEvent(iSec);
587 if ( !vTimes ) continue;
588
589 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
3cd27a08 590 Float_t time0 = fVTime0Offset[iSec]/fVTime0OffsetCounter[iSec];
591 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
467a4b25 592
3cd27a08 593 GetHistoT0(iSec,kTRUE)->Fill( time-time0,iChannel );
467a4b25 594
595
596 //Debug start
597 if ( fDebugLevel>0 ){
598 if ( !fDebugStreamer ) {
599 //debug stream
600 TDirectory *backup = gDirectory;
601 fDebugStreamer = new TTreeSRedirector("deb2.root");
602 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
603 }
604
605 Int_t row=0;
606 Int_t pad=0;
607 Int_t padc=0;
608
3cd27a08 609 Float_t q = (*GetPadQEvent(iSec)).GetMatrixArray()[iChannel];
610 Float_t rms = (*GetPadRMSEvent(iSec)).GetMatrixArray()[iChannel];
467a4b25 611
612 UInt_t channel=iChannel;
613 Int_t sector=iSec;
614
615 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
616 pad = channel-fROC->GetRowIndexes(sector)[row];
617 padc = pad-(fROC->GetNPads(sector,row)/2);
618
619 TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
620 Form("hSignalD%d.%d.%d",sector,row,pad),
621 fLastTimeBin-fFirstTimeBin,
622 fFirstTimeBin,fLastTimeBin);
623 h1->SetDirectory(0);
624
625 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
626 h1->Fill(i,fPadSignal(i));
627
628 (*fDebugStreamer) << "DataPad" <<
3cd27a08 629// "Event=" << fEvent <<
467a4b25 630 "Sector="<< sector <<
631 "Row=" << row<<
632 "Pad=" << pad <<
633 "PadC=" << padc <<
634 "PadSec="<< channel <<
3cd27a08 635 "Time0=" << time0 <<
636 "Time=" << time <<
637 "RMS=" << rms <<
638 "Sum=" << q <<
467a4b25 639 "hist.=" << h1 <<
640 "\n";
641
642 delete h1;
643 }
644 //Debug end
467a4b25 645 }
646 }
08205ed7 647}
648//_____________________________________________________________________
649Bool_t AliTPCCalibPulser::ProcessEventFast(AliTPCRawStreamFast *rawStreamFast)
650{
651 //
652 // Event Processing loop - AliTPCRawStream
653 //
654 ResetEvent();
467a4b25 655
08205ed7 656 Bool_t withInput = kFALSE;
657
658 while ( rawStreamFast->NextDDL() ){
659 while ( rawStreamFast->NextChannel() ){
660 Int_t isector = rawStreamFast->GetSector(); // current sector
661 Int_t iRow = rawStreamFast->GetRow(); // current row
662 Int_t iPad = rawStreamFast->GetPad(); // current pad
663 Int_t startTbin = (Int_t)rawStreamFast->GetStartTimeBin();
664 Int_t endTbin = (Int_t)rawStreamFast->GetEndTimeBin();
665
666 while ( rawStreamFast->NextBunch() ){
667 for (Int_t iTimeBin = startTbin; iTimeBin < endTbin; iTimeBin++){
668 Float_t signal=(Float_t)rawStreamFast->GetSignals()[iTimeBin-startTbin];
669 Update(isector,iRow,iPad,iTimeBin+1,signal);
670 withInput = kTRUE;
671 }
672 }
673 }
674 }
675 if (withInput){
676 EndEvent();
677 }
678 return withInput;
679}
680//_____________________________________________________________________
681Bool_t AliTPCCalibPulser::ProcessEventFast(AliRawReader *rawReader)
682{
683 //
684 // Event processing loop - AliRawReader
685 //
686 AliTPCRawStreamFast *rawStreamFast = new AliTPCRawStreamFast(rawReader);
687 Bool_t res=ProcessEventFast(rawStreamFast);
688 delete rawStreamFast;
689 return res;
467a4b25 690}
691//_____________________________________________________________________
3cd27a08 692Bool_t AliTPCCalibPulser::ProcessEvent(AliTPCRawStream *rawStream)
467a4b25 693{
694 //
695 // Event Processing loop - AliTPCRawStream
696 //
697
698 rawStream->SetOldRCUFormat(fOldRCUformat);
699
700 ResetEvent();
701
702 Bool_t withInput = kFALSE;
703
704 while (rawStream->Next()) {
467a4b25 705 Int_t isector = rawStream->GetSector(); // current sector
706 Int_t iRow = rawStream->GetRow(); // current row
707 Int_t iPad = rawStream->GetPad(); // current pad
708 Int_t iTimeBin = rawStream->GetTime(); // current time bin
709 Float_t signal = rawStream->GetSignal(); // current ADC signal
710
711 Update(isector,iRow,iPad,iTimeBin,signal);
712 withInput = kTRUE;
713 }
467a4b25 714 if (withInput){
715 EndEvent();
716 }
467a4b25 717 return withInput;
718}
719//_____________________________________________________________________
720Bool_t AliTPCCalibPulser::ProcessEvent(AliRawReader *rawReader)
721{
722 //
723 // Event processing loop - AliRawReader
724 //
725
726
727 AliTPCRawStream rawStream(rawReader);
728
729 rawReader->Select("TPC");
730
731 return ProcessEvent(&rawStream);
732}
733//_____________________________________________________________________
734Bool_t AliTPCCalibPulser::ProcessEvent(eventHeaderStruct *event)
735{
736 //
737 // Event processing loop - date event
738 //
739 AliRawReader *rawReader = new AliRawReaderDate((void*)event);
740 Bool_t result=ProcessEvent(rawReader);
741 delete rawReader;
742 return result;
743
744}
745//_____________________________________________________________________
3cd27a08 746TH2S* AliTPCCalibPulser::GetHisto(Int_t sector, TObjArray *arr,
467a4b25 747 Int_t nbinsY, Float_t ymin, Float_t ymax,
748 Char_t *type, Bool_t force)
749{
750 //
751 // return pointer to Q histogram
752 // if force is true create a new histogram if it doesn't exist allready
753 //
754 if ( !force || arr->UncheckedAt(sector) )
755 return (TH2S*)arr->UncheckedAt(sector);
756
757 // if we are forced and histogram doesn't yes exist create it
758 Char_t name[255], title[255];
759
760 sprintf(name,"hCalib%s%.2d",type,sector);
761 sprintf(title,"%s calibration histogram sector %.2d",type,sector);
762
763 // new histogram with Q calib information. One value for each pad!
764 TH2S* hist = new TH2S(name,title,
765 nbinsY, ymin, ymax,
766 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
767 hist->SetDirectory(0);
768 arr->AddAt(hist,sector);
769 return hist;
770}
771//_____________________________________________________________________
3cd27a08 772TH2S* AliTPCCalibPulser::GetHistoT0(Int_t sector, Bool_t force)
467a4b25 773{
774 //
775 // return pointer to T0 histogram
776 // if force is true create a new histogram if it doesn't exist allready
777 //
778 TObjArray *arr = &fHistoT0Array;
779 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
780}
781//_____________________________________________________________________
3cd27a08 782TH2S* AliTPCCalibPulser::GetHistoQ(Int_t sector, Bool_t force)
467a4b25 783{
784 //
785 // return pointer to Q histogram
786 // if force is true create a new histogram if it doesn't exist allready
787 //
788 TObjArray *arr = &fHistoQArray;
789 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
790}
791//_____________________________________________________________________
3cd27a08 792TH2S* AliTPCCalibPulser::GetHistoRMS(Int_t sector, Bool_t force)
467a4b25 793{
794 //
795 // return pointer to Q histogram
796 // if force is true create a new histogram if it doesn't exist allready
797 //
798 TObjArray *arr = &fHistoRMSArray;
799 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
800}
801//_____________________________________________________________________
3cd27a08 802TVectorF* AliTPCCalibPulser::GetPadInfoEvent(Int_t sector, TObjArray *arr, Bool_t force)
467a4b25 803{
804 //
805 // return pointer to Pad Info from 'arr' for the current event and sector
806 // if force is true create it if it doesn't exist allready
807 //
808 if ( !force || arr->UncheckedAt(sector) )
809 return (TVectorF*)arr->UncheckedAt(sector);
810
811 TVectorF *vect = new TVectorF(fROC->GetNChannels(sector));
812 arr->AddAt(vect,sector);
813 return vect;
814}
815//_____________________________________________________________________
3cd27a08 816TVectorF* AliTPCCalibPulser::GetPadTimesEvent(Int_t sector, Bool_t force)
467a4b25 817{
818 //
819 // return pointer to Pad Times Array for the current event and sector
820 // if force is true create it if it doesn't exist allready
821 //
822 TObjArray *arr = &fPadTimesArrayEvent;
823 return GetPadInfoEvent(sector,arr,force);
824}
825//_____________________________________________________________________
3cd27a08 826TVectorF* AliTPCCalibPulser::GetPadQEvent(Int_t sector, Bool_t force)
467a4b25 827{
828 //
829 // return pointer to Pad Q Array for the current event and sector
830 // if force is true create it if it doesn't exist allready
831 // for debugging purposes only
832 //
833
834 TObjArray *arr = &fPadQArrayEvent;
835 return GetPadInfoEvent(sector,arr,force);
836}
837//_____________________________________________________________________
3cd27a08 838TVectorF* AliTPCCalibPulser::GetPadRMSEvent(Int_t sector, Bool_t force)
467a4b25 839{
840 //
841 // return pointer to Pad RMS Array for the current event and sector
842 // if force is true create it if it doesn't exist allready
843 // for debugging purposes only
844 //
845 TObjArray *arr = &fPadRMSArrayEvent;
846 return GetPadInfoEvent(sector,arr,force);
847}
848//_____________________________________________________________________
3cd27a08 849TVectorF* AliTPCCalibPulser::GetPadPedestalEvent(Int_t sector, Bool_t force)
467a4b25 850{
851 //
852 // return pointer to Pad RMS Array for the current event and sector
853 // if force is true create it if it doesn't exist allready
854 // for debugging purposes only
855 //
856 TObjArray *arr = &fPadPedestalArrayEvent;
857 return GetPadInfoEvent(sector,arr,force);
858}
859//_____________________________________________________________________
3cd27a08 860AliTPCCalROC* AliTPCCalibPulser::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
467a4b25 861{
862 //
863 // return pointer to ROC Calibration
864 // if force is true create a new histogram if it doesn't exist allready
865 //
866 if ( !force || arr->UncheckedAt(sector) )
867 return (AliTPCCalROC*)arr->UncheckedAt(sector);
868
869 // if we are forced and histogram doesn't yes exist create it
870
871 // new AliTPCCalROC for T0 information. One value for each pad!
872 AliTPCCalROC *croc = new AliTPCCalROC(sector);
873 arr->AddAt(croc,sector);
874 return croc;
875}
876//_____________________________________________________________________
3cd27a08 877AliTPCCalROC* AliTPCCalibPulser::GetCalRocT0(Int_t sector, Bool_t force)
467a4b25 878{
879 //
880 // return pointer to Carge ROC Calibration
881 // if force is true create a new histogram if it doesn't exist allready
882 //
883 TObjArray *arr = &fCalRocArrayT0;
884 return GetCalRoc(sector, arr, force);
885}
886//_____________________________________________________________________
3cd27a08 887AliTPCCalROC* AliTPCCalibPulser::GetCalRocQ(Int_t sector, Bool_t force)
467a4b25 888{
889 //
890 // return pointer to T0 ROC Calibration
891 // if force is true create a new histogram if it doesn't exist allready
892 //
893 TObjArray *arr = &fCalRocArrayQ;
894 return GetCalRoc(sector, arr, force);
895}
896//_____________________________________________________________________
3cd27a08 897AliTPCCalROC* AliTPCCalibPulser::GetCalRocRMS(Int_t sector, Bool_t force)
467a4b25 898{
899 //
900 // return pointer to signal width ROC Calibration
901 // if force is true create a new histogram if it doesn't exist allready
902 //
903 TObjArray *arr = &fCalRocArrayRMS;
904 return GetCalRoc(sector, arr, force);
905}
906//_____________________________________________________________________
907AliTPCCalROC* AliTPCCalibPulser::GetCalRocOutliers(Int_t sector, Bool_t force)
908{
909 //
910 // return pointer to Outliers
911 // if force is true create a new histogram if it doesn't exist allready
912 //
913 TObjArray *arr = &fCalRocArrayOutliers;
914 return GetCalRoc(sector, arr, force);
915}
916//_____________________________________________________________________
3cd27a08 917void AliTPCCalibPulser::ResetEvent()
467a4b25 918{
919 //
920 // Reset global counters -- Should be called before each event is processed
921 //
922 fLastSector=-1;
923 fCurrentSector=-1;
924 fCurrentRow=-1;
925 fCurrentChannel=-1;
926
927 ResetPad();
928
929 fPadTimesArrayEvent.Delete();
930 fPadQArrayEvent.Delete();
931 fPadRMSArrayEvent.Delete();
932 fPadPedestalArrayEvent.Delete();
933
934 for ( Int_t i=0; i<72; ++i ){
935 fVTime0Offset[i]=0;
936 fVTime0OffsetCounter[i]=0;
937 }
938}
939//_____________________________________________________________________
3cd27a08 940void AliTPCCalibPulser::ResetPad()
467a4b25 941{
942 //
943 // Reset pad infos -- Should be called after a pad has been processed
944 //
945 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
946 fPadSignal[i] = 0;
947 fMaxTimeBin = -1;
948 fMaxPadSignal = -1;
949 fPadPedestal = -1;
950 fPadNoise = -1;
951}
952//_____________________________________________________________________
953void AliTPCCalibPulser::Merge(AliTPCCalibPulser *sig)
954{
955 //
956 // Merge reference histograms of sig to the current AliTPCCalibPulser
957 //
958
959 //merge histograms
960 for (Int_t iSec=0; iSec<72; ++iSec){
961 TH2S *hRefQmerge = sig->GetHistoQ(iSec);
962 TH2S *hRefT0merge = sig->GetHistoT0(iSec);
963 TH2S *hRefRMSmerge = sig->GetHistoRMS(iSec);
964
965
966 if ( hRefQmerge ){
967 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
968 TH2S *hRefQ = GetHistoQ(iSec);
969 if ( hRefQ ) hRefQ->Add(hRefQmerge);
970 else {
971 TH2S *hist = new TH2S(*hRefQmerge);
972 hist->SetDirectory(0);
973 fHistoQArray.AddAt(hist, iSec);
974 }
975 hRefQmerge->SetDirectory(dir);
976 }
977 if ( hRefT0merge ){
978 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
979 TH2S *hRefT0 = GetHistoT0(iSec);
980 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
981 else {
982 TH2S *hist = new TH2S(*hRefT0merge);
983 hist->SetDirectory(0);
984 fHistoT0Array.AddAt(hist, iSec);
985 }
986 hRefT0merge->SetDirectory(dir);
987 }
988 if ( hRefRMSmerge ){
989 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
990 TH2S *hRefRMS = GetHistoRMS(iSec);
991 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
992 else {
993 TH2S *hist = new TH2S(*hRefRMSmerge);
994 hist->SetDirectory(0);
995 fHistoRMSArray.AddAt(hist, iSec);
996 }
997 hRefRMSmerge->SetDirectory(dir);
998 }
999
1000 }
1001}
1002//_____________________________________________________________________
3cd27a08 1003void AliTPCCalibPulser::Analyse()
467a4b25 1004{
1005 //
1006 // Calculate calibration constants
1007 //
1008
1009 TVectorD paramQ(3);
1010 TVectorD paramT0(3);
1011 TVectorD paramRMS(3);
1012 TMatrixD dummy(3,3);
1013
1014 for (Int_t iSec=0; iSec<72; ++iSec){
1015 TH2S *hT0 = GetHistoT0(iSec);
1016 if (!hT0 ) continue;
1017
1018 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1019 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1020 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1021 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1022
1023 TH2S *hQ = GetHistoQ(iSec);
1024 TH2S *hRMS = GetHistoRMS(iSec);
1025
3cd27a08 1026 Short_t *arrayhQ = hQ->GetArray();
1027 Short_t *arrayhT0 = hT0->GetArray();
1028 Short_t *arrayhRMS = hRMS->GetArray();
467a4b25 1029
1030 UInt_t nChannels = fROC->GetNChannels(iSec);
1031
1032 //debug
1033 Int_t row=0;
1034 Int_t pad=0;
1035 Int_t padc=0;
1036 //! debug
1037
1038 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1039
467a4b25 1040 Float_t cogTime0 = -1000;
1041 Float_t cogQ = -1000;
1042 Float_t cogRMS = -1000;
1043 Float_t cogOut = 0;
1044
467a4b25 1045 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1046 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1047 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
467a4b25 1048/*
3cd27a08 1049 AliMathBase::FitGaus(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&paramQ,&dummy);
1050 AliMathBase::FitGaus(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&paramT0,&dummy);
1051 AliMathBase::FitGaus(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&paramRMS,&dummy);
467a4b25 1052 cogQ = paramQ[1];
1053 cogTime0 = paramT0[1];
1054 cogRMS = paramRMS[1];
1055*/
3cd27a08 1056 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ);
1057 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0);
1058 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS);
467a4b25 1059
467a4b25 1060 /*
1061 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
1062 cogOut = 1;
1063 cogTime0 = 0;
1064 cogQ = 0;
1065 cogRMS = 0;
1066 }
1067*/
1068 rocQ->SetValue(iChannel, cogQ*cogQ);
1069 rocT0->SetValue(iChannel, cogTime0);
1070 rocRMS->SetValue(iChannel, cogRMS);
1071 rocOut->SetValue(iChannel, cogOut);
1072
1073
1074 //debug
1075 if ( fDebugLevel > 0 ){
1076 if ( !fDebugStreamer ) {
1077 //debug stream
1078 TDirectory *backup = gDirectory;
1079 fDebugStreamer = new TTreeSRedirector("deb2.root");
1080 if ( backup ) backup->cd(); //we don't want to be cd'd to the debug streamer
1081 }
1082
1083 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
1084 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
1085 padc = pad-(fROC->GetNPads(iSec,row)/2);
1086
1087 (*fDebugStreamer) << "DataEnd" <<
1088 "Sector=" << iSec <<
1089 "Pad=" << pad <<
1090 "PadC=" << padc <<
1091 "Row=" << row <<
1092 "PadSec=" << iChannel <<
1093 "Q=" << cogQ <<
1094 "T0=" << cogTime0 <<
1095 "RMS=" << cogRMS <<
1096 "\n";
1097 }
1098 //! debug
467a4b25 1099 }
467a4b25 1100 }
1101 delete fDebugStreamer;
1102 fDebugStreamer = 0x0;
1103}
1104//_____________________________________________________________________
1105void AliTPCCalibPulser::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t append)
1106{
1107 //
1108 // Write class to file
1109 //
1110
1111 TString sDir(dir);
1112 TString option;
1113
1114 if ( append )
1115 option = "update";
1116 else
1117 option = "recreate";
1118
1119 TDirectory *backup = gDirectory;
1120 TFile f(filename,option.Data());
1121 f.cd();
1122 if ( !sDir.IsNull() ){
1123 f.mkdir(sDir.Data());
1124 f.cd(sDir);
1125 }
1126 this->Write();
1127 f.Close();
1128
1129 if ( backup ) backup->cd();
1130}
1131//_____________________________________________________________________
1132//_________________________ Test Functions ___________________________
1133//_____________________________________________________________________
1134TObjArray* AliTPCCalibPulser::TestBinning()
1135{
1136 //
1137 // Function to test the binning of the reference histograms
1138 // type: T0, Q or RMS
1139 // mode: 0 - number of filled bins per channel
1140 // 1 - number of empty bins between filled bins in one ROC
1141 // returns TObjArray with the test histograms type*2+mode:
1142 // position 0 = T0,0 ; 1 = T0,1 ; 2 = Q,0 ...
1143
1144
1145 TObjArray *histArray = new TObjArray(6);
1146 const Char_t *type[] = {"T0","Q","RMS"};
1147 Int_t fNbins[3] = {fNbinsT0,fNbinsQ,fNbinsRMS};
1148
1149 for (Int_t itype = 0; itype<3; ++itype){
1150 for (Int_t imode=0; imode<2; ++imode){
1151 Int_t icount = itype*2+imode;
1152 histArray->AddAt(new TH1F(Form("hTestBinning%s%d",type[itype],imode),
1153 Form("Test Binning of '%s', mode - %d",type[itype],imode),
1154 72,0,72),
1155 icount);
1156 }
1157 }
1158
1159
1160 TH2S *hRef=0x0;
1161 Short_t *array=0x0;
1162 for (Int_t itype = 0; itype<3; ++itype){
1163 for (Int_t iSec=0; iSec<72; ++iSec){
1164 if ( itype == 0 ) hRef = GetHistoT0(iSec);
1165 if ( itype == 1 ) hRef = GetHistoQ(iSec);
1166 if ( itype == 2 ) hRef = GetHistoRMS(iSec);
1167 if ( hRef == 0x0 ) continue;
1168 array = (hRef->GetArray());
1169 UInt_t nChannels = fROC->GetNChannels(iSec);
1170
1171 Int_t nempty=0;
1172 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1173 Int_t nfilled=0;
1174 Int_t offset = (fNbins[itype]+2)*(iChannel+1)+1;
1175 Int_t c1 = 0;
1176 Int_t c2 = 0;
1177 for (Int_t iBin=0; iBin<fNbins[itype]; ++iBin){
1178 if ( array[offset+iBin]>0 ) {
1179 nfilled++;
1180 if ( c1 && c2 ) nempty++;
1181 else c1 = 1;
1182 }
1183 else if ( c1 ) c2 = 1;
1184
1185 }
1186 ((TH1F*)histArray->At(itype*2))->Fill(nfilled);
1187 }
1188 ((TH1F*)histArray->At(itype*2+1))->Fill(iSec,nempty);
1189 }
1190 }
1191 return histArray;
1192}