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