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