]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TPC/AliTPCCalibCE.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing Coverity 10889 (B.Hippolyte)
[u/mrichter/AliRoot.git] / TPC / AliTPCCalibCE.cxx
CommitLineData
75d8233f 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
75d8233f 16/* $Id$ */
17
3cd27a08 18////////////////////////////////////////////////////////////////////////////////////////
19// //
20// Implementation of the TPC Central Electrode calibration //
21// //
22// Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
23// //
24////////////////////////////////////////////////////////////////////////////////////////
7fb602b1 25//
26//
27// *************************************************************************************
28// * Class Description *
29// *************************************************************************************
30//
31/* BEGIN_HTML
32 <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
33 using laser runs.</h4>
34
35 The information retrieved is
36 <ul style="list-style-type: square;">
37 <li>Time arrival from the CE</li>
38 <li>Signal width</li>
39 <li>Signal sum</li>
40 </ul>
41
42<h4>Overview:</h4>
43 <ol style="list-style-type: upper-roman;">
44 <li><a href="#working">Working principle</a></li>
45 <li><a href="#user">User interface for filling data</a></li>
46 <li><a href="#info">Stored information</a></li>
47 </ol>
48
49 <h3><a name="working">I. Working principle</a></h3>
50
51 <h4>Raw laser data is processed by calling one of the ProcessEvent(...) functions
52 (see below). These in the end call the Update(...) function.</h4>
53
54 <ul style="list-style-type: square;">
55 <li>the Update(...) function:<br />
56 In this function the array fPadSignal is filled with the adc signals between the specified range
57 fFirstTimeBin and fLastTimeBin for the current pad.
58 before going to the next pad the ProcessPad() function is called, which analyses the data for one pad
59 stored in fPadSignal.
60 </li>
61 <ul style="list-style-type: square;">
62 <li>the ProcessPad() function:</li>
63 <ol style="list-style-type: decimal;">
64 <li>Find Pedestal and Noise information</li>
65 <ul style="list-style-type: square;">
66 <li>use database information which has to be set by calling<br />
67 SetPedestalDatabase(AliTPCCalPad *pedestalTPC, AliTPCCalPad *padNoiseTPC)</li>
68 <li>if no information from the pedestal data base
69 is available the informaion is calculated on the fly
70 ( see FindPedestal() function )</li>
71 </ul>
72 <li>Find local maxima of the pad signal</li>
73 <ul style="list-style-type: square;">
74 <li>maxima arise from the laser tracks, the CE and also periodic postpeaks after the CE signal have
75 have been observed ( see FindLocalMaxima(...) )</li>
76 </ul>
77 <li>Find the CE signal information</li>
78 <ul style="list-style-type: square;">
79 <li>to find the position of the CE signal the Tmean information from the previos event is used
80 as the CE signal the local maximum closest to this Tmean is identified</li>
81 <li>calculate mean = T0, RMS = signal width and Q sum in a range of -4+7 timebins around Q max position
82 the Q sum is scaled by pad area (see FindPulserSignal(...) function)</li>
83 </ul>
84 <li>Fill a temprary array for the T0 information (GetPadTimesEvent(fCurrentSector,kTRUE)) (why see below)</li>
85 <li>Fill the Q sum and RMS values in the histograms (GetHisto[RMS,Q](ROC,kTRUE))</li>
86 </ol>
87 </ul>
88 </ul>
89
90 <h4>At the end of each event the EndEvent() function is called</h4>
91
92 <ul style="list-style-type: square;">
93 <li>the EndEvent() function:</li>
94 <ul style="list-style-type: square;">
95 <li>calculate the mean T0 for side A and side C. Fill T0 histogram with Time0-<Time0 for side[A,C]>
96 This is done to overcome syncronisation problems between the trigger and the fec clock.</li>
97 <li>calculate Mean T for each ROC using the COG aroud the median of the LocalMaxima distribution in one sector</li>
98 <li>calculate Mean Q</li>
99 <li>calculate Global fit parameters for Pol1 and Pol2 fits</li>
100 </ul>
101 </ul>
102
103 <h4>After accumulating the desired statistics the Analyse() function has to be called.</h4>
104 <ul style="list-style-type: square;">
105 <li>the Analyse() function:</li>
106 <ul style="list-style-type: square;">
107 <li>calculate the mean values of T0, RMS, Q for each pad, using
108 the AliMathBase::GetCOG(...) function</li>
109 <li>fill the calibration storage classes (AliTPCCalROC) for each ROC</li>
110 (The calibration information is stored in the TObjArrays fCalRocArrayT0, fCalRocArrayRMS and fCalRocArrayQ</li>
111 </ul>
112 </ul>
113
114 <h3><a name="user">II. User interface for filling data</a></h3>
115
116 <h4>To Fill information one of the following functions can be used:</h4>
117
118 <ul style="list-style-type: none;">
119 <li> Bool_t ProcessEvent(eventHeaderStruct *event);</li>
120 <ul style="list-style-type: square;">
121 <li>process Date event</li>
122 <li>use AliTPCRawReaderDate and call ProcessEvent(AliRawReader *rawReader)</li>
123 </ul>
124 <br />
125
126 <li> Bool_t ProcessEvent(AliRawReader *rawReader);</li>
127 <ul style="list-style-type: square;">
128 <li>process AliRawReader event</li>
129 <li>use AliTPCRawStream to loop over data and call ProcessEvent(AliTPCRawStream *rawStream)</li>
130 </ul>
131 <br />
132
133 <li> Bool_t ProcessEvent(AliTPCRawStream *rawStream);</li>
134 <ul style="list-style-type: square;">
135 <li>process event from AliTPCRawStream</li>
136 <li>call Update function for signal filling</li>
137 </ul>
138 <br />
139
140 <li> Int_t Update(const Int_t isector, const Int_t iRow, const Int_t
141 iPad, const Int_t iTimeBin, const Float_t signal);</li>
142 <ul style="list-style-type: square;">
143 <li>directly fill signal information (sector, row, pad, time bin, pad)
144 to the reference histograms</li>
145 </ul>
146 </ul>
147
148 <h4>It is also possible to merge two independently taken calibrations using the function</h4>
149
150 <ul style="list-style-type: none;">
151 <li> void Merge(AliTPCCalibSignal *sig)</li>
152 <ul style="list-style-type: square;">
153 <li>copy histograms in 'sig' if they do not exist in this instance</li>
154 <li>Add histograms in 'sig' to the histograms in this instance if the allready exist</li>
155 <li>After merging call Analyse again!</li>
156 </ul>
157 </ul>
158
159
160 <h4>example: filling data using root raw data:</h4>
161 <pre>
162 void fillCE(Char_t *filename)
163 {
164 rawReader = new AliRawReaderRoot(fileName);
165 if ( !rawReader ) return;
166 AliTPCCalibCE *calib = new AliTPCCalibCE;
167 while (rawReader->NextEvent()){
168 calib->ProcessEvent(rawReader);
169 }
170 calib->Analyse();
171 calib->DumpToFile("CEData.root");
172 delete rawReader;
173 delete calib;
174 }
175 </pre>
176
177 <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
178
179 <h4><a name="info:stored">III.1 Stored information</a></h4>
180 <ul style="list-style-type: none;">
181 <li>Histograms:</li>
182 <ul style="list-style-type: none;">
183 <li>For each ROC three TH2S histos 'Reference Histograms' (ROC channel vs. [Time0, signal width, Q sum])
184 is created when it is filled for the first time (GetHisto[T0,RMS,Q](ROC,kTRUE)). The histos are
185 stored in the TObjArrays fHistoT0Array, fHistoRMSArray and fHistoQArray.</li>
186 </ul>
187 <br />
188
189 <li>Calibration Data:</li>
190 <ul style="list-style-type: none;">
191 <li>For each ROC three types of calibration data (AliTPCCalROC) is stored: for the mean arrival Time,
192 the signal width and the signal Sum. The AliTPCCalROC objects are stored in the TObjArrays
193 fCalRocArrayT0, fCalRocArrayRMS , fCalRocArrayQ. The object for each roc is created the first time it
194 is accessed (GetCalRoc[T0,RMS,Q](ROC,kTRUE));</li>
195 </ul>
196 <br />
197
198 <li>For each event the following information is stored:</li>
199
200 <ul style="list-style-type: square;">
201 <li>event time ( TVectorD fVEventTime )</li>
202 <li>event id ( TVectorD fVEventNumber )</li>
203 <br />
204 <li>mean arrival time for each ROC ( TObjArray fTMeanArrayEvent )</li>
205 <li>mean Q for each ROC ( TObjArray fQMeanArrayEvent )</li>
206 <li>parameters of a plane fit for each ROC ( TObjArray fParamArrayEventPol1 )</li>
207 <li>parameters of a 2D parabola fit for each ROC ( TObjArray fParamArrayEventPol2 )</li>
208 </ul>
209 </ul>
210
211 <h4><a name="info:retrieve">III.2 Retrieving information</a></h4>
212 <ul style="list-style-type: none;">
213 <li>Accessing the 'Reference Histograms' (Time0, signal width and Q sum information pad by pad):</li>
214 <ul style="list-style-type: square;">
215 <li>TH2F *GetHistoT0(Int_t sector);</li>
216 <li>TH2F *GetHistoRMS(Int_t sector);</li>
217 <li>TH2F *GetHistoQ(Int_t sector);</li>
218 </ul>
219 <br />
220
221 <li>Accessing the calibration storage objects:</li>
222 <ul style="list-style-type: square;">
223 <li>AliTPCCalROC *GetCalRocT0(Int_t sector); // for the Time0 values</li>
224 <li>AliTPCCalROC *GetCalRocRMS(Int_t sector); // for the signal width values</li>
225 <li>AliTPCCalROC *GetCalRocQ(Int_t sector); // for the Q sum values</li>
226 </ul>
227 <br />
228
229 <li>Accessin the event by event information:</li>
230 <ul style="list-style-type: square;">
231 <li>The event by event information can be displayed using the</li>
232 <li>MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)</li>
233 <li>which creates a graph from the specified variables</li>
234 </ul>
235 </ul>
236
237 <h4>example for visualisation:</h4>
238 <pre>
239 //if the file "CEData.root" was created using the above example one could do the following:
240 TFile fileCE("CEData.root")
241 AliTPCCalibCE *ce = (AliTPCCalibCE*)fileCE->Get("AliTPCCalibCE");
242 ce->GetCalRocT0(0)->Draw("colz");
243 ce->GetCalRocRMS(0)->Draw("colz");
244
245 //or use the AliTPCCalPad functionality:
246 AliTPCCalPad padT0(ped->GetCalPadT0());
247 AliTPCCalPad padSigWidth(ped->GetCalPadRMS());
248 padT0->MakeHisto2D()->Draw("colz"); //Draw A-Side Time0 Information
249 padSigWidth->MakeHisto2D()->Draw("colz"); //Draw A-Side signal width Information
250
251 //display event by event information:
252 //Draw mean arrival time as a function of the event time for oroc sector A00
253 ce->MakeGraphTimeCE(36, 0, 2)->Draw("alp");
254 //Draw first derivative in local x from a plane fit as a function of the event time for oroc sector A00
255 ce->MakeGraphTimeCE(36, 0, 0, 1)->Draw("alp");
256 </pre>
257END_HTML */
258//////////////////////////////////////////////////////////////////////////////////////
259
260
3cd27a08 261//Root includes
262#include <TObjArray.h>
78f17711 263#include <TH1.h>
3cd27a08 264#include <TH1F.h>
265#include <TH2S.h>
78f17711 266#include <TF1.h>
3cd27a08 267#include <TString.h>
268#include <TVectorF.h>
269#include <TVectorD.h>
78f17711 270#include <TVector3.h>
3cd27a08 271#include <TMatrixD.h>
272#include <TMath.h>
273#include <TGraph.h>
78f17711 274#include <TGraphErrors.h>
3cd27a08 275#include <TString.h>
ac940b58 276#include <TMap.h>
3cd27a08 277#include <TDirectory.h>
278#include <TSystem.h>
279#include <TFile.h>
7442bceb 280#include <TCollection.h>
78f17711 281#include <TTimeStamp.h>
282#include <TList.h>
283#include <TKey.h>
3cd27a08 284
285//AliRoot includes
4c6d06dc 286#include "AliLog.h"
3cd27a08 287#include "AliRawReader.h"
288#include "AliRawReaderRoot.h"
289#include "AliRawReaderDate.h"
290#include "AliRawEventHeaderBase.h"
291#include "AliTPCRawStream.h"
3cd27a08 292#include "AliTPCCalROC.h"
293#include "AliTPCCalPad.h"
294#include "AliTPCROC.h"
295#include "AliTPCParam.h"
296#include "AliTPCCalibCE.h"
297#include "AliMathBase.h"
78f17711 298#include "AliTPCTransform.h"
299#include "AliTPCLaserTrack.h"
3cd27a08 300#include "TTreeStream.h"
301
78f17711 302#include "AliCDBManager.h"
303#include "AliCDBEntry.h"
3cd27a08 304//date
305#include "event.h"
306ClassImp(AliTPCCalibCE)
307
308
7fb602b1 309AliTPCCalibCE::AliTPCCalibCE() :
880c3382 310 AliTPCCalibRawBase(),
311 fNbinsT0(200),
312 fXminT0(-5),
313 fXmaxT0(5),
314 fNbinsQ(200),
315 fXminQ(1),
316 fXmaxQ(40),
317 fNbinsRMS(100),
318 fXminRMS(0.1),
319 fXmaxRMS(5.1),
320 fPeakDetMinus(2),
321 fPeakDetPlus(3),
322 fPeakIntMinus(2),
323 fPeakIntPlus(2),
324 fNoiseThresholdMax(5.),
325 fNoiseThresholdSum(8.),
326 fIsZeroSuppressed(kFALSE),
327 fLastSector(-1),
328 fSecRejectRatio(.4),
329 fParam(new AliTPCParam),
330 fPedestalTPC(0x0),
331 fPadNoiseTPC(0x0),
332 fPedestalROC(0x0),
333 fPadNoiseROC(0x0),
334 fCalRocArrayT0(72),
335 fCalRocArrayT0Err(72),
336 fCalRocArrayQ(72),
337 fCalRocArrayRMS(72),
338 fCalRocArrayOutliers(72),
339 fHistoQArray(72),
340 fHistoT0Array(72),
341 fHistoRMSArray(72),
342 fMeanT0rms(0),
343 fMeanQrms(0),
344 fMeanRMSrms(0),
345 fHistoTmean(72),
346 fParamArrayEventPol1(72),
347 fParamArrayEventPol2(72),
348 fTMeanArrayEvent(72),
349 fQMeanArrayEvent(72),
350 fVEventTime(1000),
351 fVEventNumber(1000),
352 fVTime0SideA(1000),
353 fVTime0SideC(1000),
880c3382 354 fEventId(-1),
c3066940 355 fOldRunNumber(0),
880c3382 356 fPadTimesArrayEvent(72),
357 fPadQArrayEvent(72),
358 fPadRMSArrayEvent(72),
359 fPadPedestalArrayEvent(72),
360 fCurrentChannel(-1),
361 fCurrentSector(-1),
362 fCurrentRow(-1),
363 fMaxPadSignal(-1),
364 fMaxTimeBin(-1),
2963bcbf 365// fPadSignal(1024),
880c3382 366 fPadPedestal(0),
367 fPadNoise(0),
368 fVTime0Offset(72),
369 fVTime0OffsetCounter(72),
370 fVMeanQ(72),
371 fVMeanQCounter(72),
78f17711 372 fCurrentCETimeRef(0),
373 fProcessOld(kTRUE),
374 fProcessNew(kFALSE),
375 fAnalyseNew(kTRUE),
376 fHnDrift(0x0),
377 fArrHnDrift(100),
378 fTimeBursts(100),
379 fArrFitGraphs(0x0)
75d8233f 380{
880c3382 381 //
382 // AliTPCSignal default constructor
383 //
384 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
385 fFirstTimeBin=650;
78f17711 386 fLastTimeBin=1030;
880c3382 387 fParam->Update();
2963bcbf 388 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
78f17711 389 for (Int_t i=0;i<5;++i){
390 fPeaks[i]=0;
391 fPeakWidths[i]=0;
392 }
393 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
75d8233f 394}
395//_____________________________________________________________________
396AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
880c3382 397 AliTPCCalibRawBase(sig),
398 fNbinsT0(sig.fNbinsT0),
399 fXminT0(sig.fXminT0),
400 fXmaxT0(sig.fXmaxT0),
401 fNbinsQ(sig.fNbinsQ),
402 fXminQ(sig.fXminQ),
403 fXmaxQ(sig.fXmaxQ),
404 fNbinsRMS(sig.fNbinsRMS),
405 fXminRMS(sig.fXminRMS),
406 fXmaxRMS(sig.fXmaxRMS),
407 fPeakDetMinus(sig.fPeakDetMinus),
408 fPeakDetPlus(sig.fPeakDetPlus),
409 fPeakIntMinus(sig.fPeakIntMinus),
410 fPeakIntPlus(sig.fPeakIntPlus),
411 fNoiseThresholdMax(sig.fNoiseThresholdMax),
412 fNoiseThresholdSum(sig.fNoiseThresholdSum),
413 fIsZeroSuppressed(sig.fIsZeroSuppressed),
414 fLastSector(-1),
415 fSecRejectRatio(.4),
416 fParam(new AliTPCParam),
417 fPedestalTPC(0x0),
418 fPadNoiseTPC(0x0),
419 fPedestalROC(0x0),
420 fPadNoiseROC(0x0),
421 fCalRocArrayT0(72),
422 fCalRocArrayT0Err(72),
423 fCalRocArrayQ(72),
424 fCalRocArrayRMS(72),
425 fCalRocArrayOutliers(72),
426 fHistoQArray(72),
427 fHistoT0Array(72),
428 fHistoRMSArray(72),
429 fMeanT0rms(sig.fMeanT0rms),
430 fMeanQrms(sig.fMeanQrms),
431 fMeanRMSrms(sig.fMeanRMSrms),
432 fHistoTmean(72),
433 fParamArrayEventPol1(72),
434 fParamArrayEventPol2(72),
435 fTMeanArrayEvent(72),
436 fQMeanArrayEvent(72),
437 fVEventTime(sig.fVEventTime),
438 fVEventNumber(sig.fVEventNumber),
439 fVTime0SideA(sig.fVTime0SideA),
440 fVTime0SideC(sig.fVTime0SideC),
880c3382 441 fEventId(-1),
c3066940 442 fOldRunNumber(0),
880c3382 443 fPadTimesArrayEvent(72),
444 fPadQArrayEvent(72),
445 fPadRMSArrayEvent(72),
446 fPadPedestalArrayEvent(72),
447 fCurrentChannel(-1),
448 fCurrentSector(-1),
449 fCurrentRow(-1),
450 fMaxPadSignal(-1),
451 fMaxTimeBin(-1),
2963bcbf 452// fPadSignal(1024),
880c3382 453 fPadPedestal(0),
454 fPadNoise(0),
455 fVTime0Offset(72),
456 fVTime0OffsetCounter(72),
457 fVMeanQ(72),
458 fVMeanQCounter(72),
78f17711 459 fCurrentCETimeRef(0),
460 fProcessOld(sig.fProcessOld),
461 fProcessNew(sig.fProcessNew),
462 fAnalyseNew(sig.fAnalyseNew),
463 fHnDrift(0x0),
464 fArrHnDrift(100),
465 fTimeBursts(100),
466 fArrFitGraphs(0x0)
75d8233f 467{
ac940b58 468 //
880c3382 469 // AliTPCSignal copy constructor
ac940b58 470 //
2963bcbf 471 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
472
ac940b58 473 for (Int_t iSec = 0; iSec < 72; ++iSec){
474 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
475 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
476 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
477 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
75d8233f 478
ac940b58 479 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
480 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
481 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
75d8233f 482
ac940b58 483 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
484 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
485 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
486 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
75d8233f 487
ac940b58 488 if ( hQ != 0x0 ){
ac940b58 489 TH2S *hNew = new TH2S(*hQ);
490 hNew->SetDirectory(0);
491 fHistoQArray.AddAt(hNew,iSec);
ac940b58 492 }
493 if ( hT0 != 0x0 ){
ac940b58 494 TH2S *hNew = new TH2S(*hT0);
495 hNew->SetDirectory(0);
496 fHistoT0Array.AddAt(hNew,iSec);
ac940b58 497 }
498 if ( hRMS != 0x0 ){
ac940b58 499 TH2S *hNew = new TH2S(*hRMS);
500 hNew->SetDirectory(0);
501 fHistoRMSArray.AddAt(hNew,iSec);
75d8233f 502 }
ac940b58 503 }
75d8233f 504
880c3382 505 //copy fit parameters event by event
ac940b58 506 TObjArray *arr=0x0;
507 for (Int_t iSec=0; iSec<72; ++iSec){
508 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
509 if ( arr ){
510 TObjArray *arrEvents = new TObjArray(arr->GetSize());
511 fParamArrayEventPol1.AddAt(arrEvents, iSec);
512 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
513 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
514 arrEvents->AddAt(new TVectorD(*vec),iEvent);
515 }
7fb602b1 516
ac940b58 517 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
518 if ( arr ){
519 TObjArray *arrEvents = new TObjArray(arr->GetSize());
520 fParamArrayEventPol2.AddAt(arrEvents, iSec);
521 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
522 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
523 arrEvents->AddAt(new TVectorD(*vec),iEvent);
7fb602b1 524 }
525
ac940b58 526 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
527 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
528 if ( vMeanTime )
529 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
530 if ( vMeanQ )
531 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
532 }
533
7fb602b1 534
ac940b58 535 fVEventTime.ResizeTo(sig.fVEventTime);
536 fVEventNumber.ResizeTo(sig.fVEventNumber);
537 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
538 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
7fb602b1 539
ac940b58 540 fParam->Update();
78f17711 541
542 for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
543 TObject *o=sig.fArrHnDrift.UncheckedAt(i);
544 if (o){
545 TObject *newo=o->Clone("fHnDrift");
546 fArrHnDrift.AddAt(newo,i);
547 if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
548 }
549 }
550
551 for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
552 fTimeBursts[i]=sig.fTimeBursts[i];
553 }
554
555 for (Int_t i=0;i<5;++i){
556 fPeaks[i]=sig.fPeaks[i];
557 fPeakWidths[i]=sig.fPeakWidths[i];
558 }
559 if (sig.fArrFitGraphs) {
560 fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
561 fArrFitGraphs->SetOwner();
562 }
563
564 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
565
75d8233f 566}
ac940b58 567//_____________________________________________________________________
568AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
880c3382 569 AliTPCCalibRawBase(),
ac940b58 570 fNbinsT0(200),
571 fXminT0(-5),
572 fXmaxT0(5),
573 fNbinsQ(200),
574 fXminQ(1),
575 fXmaxQ(40),
576 fNbinsRMS(100),
577 fXminRMS(0.1),
578 fXmaxRMS(5.1),
880c3382 579 fPeakDetMinus(2),
580 fPeakDetPlus(3),
581 fPeakIntMinus(2),
582 fPeakIntPlus(2),
ac940b58 583 fNoiseThresholdMax(5.),
584 fNoiseThresholdSum(8.),
585 fIsZeroSuppressed(kFALSE),
586 fLastSector(-1),
587 fSecRejectRatio(.4),
ac940b58 588 fParam(new AliTPCParam),
589 fPedestalTPC(0x0),
590 fPadNoiseTPC(0x0),
591 fPedestalROC(0x0),
592 fPadNoiseROC(0x0),
593 fCalRocArrayT0(72),
594 fCalRocArrayT0Err(72),
595 fCalRocArrayQ(72),
596 fCalRocArrayRMS(72),
597 fCalRocArrayOutliers(72),
598 fHistoQArray(72),
599 fHistoT0Array(72),
600 fHistoRMSArray(72),
601 fMeanT0rms(0),
602 fMeanQrms(0),
603 fMeanRMSrms(0),
604 fHistoTmean(72),
605 fParamArrayEventPol1(72),
606 fParamArrayEventPol2(72),
607 fTMeanArrayEvent(72),
608 fQMeanArrayEvent(72),
880c3382 609 fVEventTime(1000),
610 fVEventNumber(1000),
611 fVTime0SideA(1000),
612 fVTime0SideC(1000),
ac940b58 613 fEventId(-1),
c3066940 614 fOldRunNumber(0),
ac940b58 615 fPadTimesArrayEvent(72),
616 fPadQArrayEvent(72),
617 fPadRMSArrayEvent(72),
618 fPadPedestalArrayEvent(72),
619 fCurrentChannel(-1),
620 fCurrentSector(-1),
621 fCurrentRow(-1),
622 fMaxPadSignal(-1),
623 fMaxTimeBin(-1),
2963bcbf 624// fPadSignal(1024),
ac940b58 625 fPadPedestal(0),
626 fPadNoise(0),
627 fVTime0Offset(72),
628 fVTime0OffsetCounter(72),
629 fVMeanQ(72),
630 fVMeanQCounter(72),
78f17711 631 fCurrentCETimeRef(0),
632 fProcessOld(kTRUE),
633 fProcessNew(kFALSE),
634 fAnalyseNew(kTRUE),
635 fHnDrift(0x0),
636 fArrHnDrift(100),
637 fTimeBursts(100),
638 fArrFitGraphs(0x0)
ac940b58 639{
640 //
641 // constructor which uses a tmap as input to set some specific parameters
642 //
880c3382 643 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
644 fFirstTimeBin=650;
78f17711 645 fLastTimeBin=1030;
ac940b58 646 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
647 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
648 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
649 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
650 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
651 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
652 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
653 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
654 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
655 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
656 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
880c3382 657 if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
658 if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
659 if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
660 if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
ac940b58 661 if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
662 if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
663 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
880c3382 664 if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
ac940b58 665 if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
666
78f17711 667 if (config->GetValue("ProcessOld")) fProcessOld = (Bool_t)((TObjString*)config->GetValue("ProcessOld"))->GetString().Atoi();
668 if (config->GetValue("ProcessNew")) fProcessNew = (Bool_t)((TObjString*)config->GetValue("ProcessNew"))->GetString().Atoi();
669 if (config->GetValue("AnalyseNew")) fAnalyseNew = (Bool_t)((TObjString*)config->GetValue("AnalyseNew"))->GetString().Atoi();
670
2963bcbf 671 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
78f17711 672 for (Int_t i=0;i<5;++i){
673 fPeaks[i]=0;
674 fPeakWidths[i]=0;
675 }
2963bcbf 676
ac940b58 677 fParam->Update();
78f17711 678 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
ac940b58 679}
680
75d8233f 681//_____________________________________________________________________
682AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
683{
684 //
685 // assignment operator
686 //
687 if (&source == this) return *this;
688 new (this) AliTPCCalibCE(source);
689
690 return *this;
691}
692//_____________________________________________________________________
693AliTPCCalibCE::~AliTPCCalibCE()
694{
78f17711 695 //
696 // destructor
697 //
698
699 fCalRocArrayT0.Delete();
700 fCalRocArrayT0Err.Delete();
701 fCalRocArrayQ.Delete();
702 fCalRocArrayRMS.Delete();
703 fCalRocArrayOutliers.Delete();
704
705 fHistoQArray.Delete();
706 fHistoT0Array.Delete();
707 fHistoRMSArray.Delete();
708
709 fHistoTmean.Delete();
710
711 fParamArrayEventPol1.Delete();
712 fParamArrayEventPol2.Delete();
713 fTMeanArrayEvent.Delete();
714 fQMeanArrayEvent.Delete();
715
716 fPadTimesArrayEvent.Delete();
717 fPadQArrayEvent.Delete();
718 fPadRMSArrayEvent.Delete();
719 fPadPedestalArrayEvent.Delete();
720
721 fArrHnDrift.SetOwner();
722 fArrHnDrift.Delete();
723
724 if (fArrFitGraphs){
725 fArrFitGraphs->SetOwner();
726 delete fArrFitGraphs;
727 }
75d8233f 728}
729//_____________________________________________________________________
7fb602b1 730Int_t AliTPCCalibCE::Update(const Int_t icsector,
75d8233f 731 const Int_t icRow,
732 const Int_t icPad,
733 const Int_t icTimeBin,
734 const Float_t csignal)
735{
ac940b58 736 //
737 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
738 // no extra analysis necessary. Assumes knowledge of the signal shape!
739 // assumes that it is looped over consecutive time bins of one pad
740 //
4c6d06dc 741
78f17711 742 if (!fProcessOld) return 0;
ac940b58 743 //temp
4c6d06dc 744
b401648b 745 if (icRow<0) return 0;
746 if (icPad<0) return 0;
747 if (icTimeBin<0) return 0;
ac940b58 748 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
75d8233f 749
ac940b58 750 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
75d8233f 751
ac940b58 752 //init first pad and sector in this event
753 if ( fCurrentChannel == -1 ) {
754 fLastSector=-1;
755 fCurrentChannel = iChannel;
756 fCurrentSector = icsector;
757 fCurrentRow = icRow;
758 }
75d8233f 759
ac940b58 760 //process last pad if we change to a new one
761 if ( iChannel != fCurrentChannel ){
762 ProcessPad();
763 fLastSector=fCurrentSector;
764 fCurrentChannel = iChannel;
765 fCurrentSector = icsector;
766 fCurrentRow = icRow;
767 }
75d8233f 768
ac940b58 769 //fill signals for current pad
2963bcbf 770 fPadSignal[icTimeBin]=csignal;
ac940b58 771 if ( csignal > fMaxPadSignal ){
772 fMaxPadSignal = csignal;
773 fMaxTimeBin = icTimeBin;
774 }
775 return 0;
75d8233f 776}
78f17711 777
778//_____________________________________________________________________
779void AliTPCCalibCE::ProcessBunch(const Int_t sector, const Int_t row, const Int_t pad,
780 const Int_t length, const UInt_t startTimeBin, const UShort_t* signal)
781{
782 //
783 // new filling method to fill the THnSparse histogram
784 //
785
786 //only in new processing mode
787 if (!fProcessNew) return;
788 //don't use the IROCs
789 if (sector<36) return;
790 //only bunches with reasonable length
791 if (length<3||length>10) return;
792
793 UShort_t timeBin = (UShort_t)startTimeBin;
794 //skip first laser layer
795 if (timeBin<200) return;
796 Double_t timeBurst=SetBurstHnDrift();
797
798 //after 1 event setup peak ranges
799 if (fNevents==1&&fPeaks[4]==0) {
800 // set time range
801 fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
802 FindLaserLayers();
803 // set time range
804 fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
805 }
806
807 // After the first event only fill every 5th bin in a row with the CE information
808 Int_t padFill=pad;
809 if (fPeaks[4]>100&&TMath::Abs((Short_t)fPeaks[4]-(Short_t)timeBin)<(Short_t)fPeakWidths[4]){
810 Int_t mod=5;
811 Int_t n=pad/mod;
812 padFill=mod*n+mod/2;
813 }
814
815 //noise removal
816 if (!IsPeakInRange(timeBin+length/2)) return;
817
818 Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
819 (Double_t)padFill,(Double_t)timeBin,timeBurst};
820
821 for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
822 Float_t sig=(Float_t)signal[iTimeBin];
823 x[3]=timeBin;
824 fHnDrift->Fill(x,sig);
825 --timeBin;
826 }
827}
828//_____________________________________________________________________
829void AliTPCCalibCE::FindLaserLayers()
830{
831 //
832 // Find the laser layer positoins
833 //
834
835
836 //find CE signal position and width
837 TH1D *hproj=fHnDrift->Projection(3);
838 hproj->GetXaxis()->SetRangeUser(700,1030);
839 Int_t maxbin=hproj->GetMaximumBin();
840 Double_t binc=hproj->GetBinCenter(maxbin);
841 hproj->GetXaxis()->SetRangeUser(binc-10,binc+10);
842
843 fPeaks[4]=(UShort_t)TMath::Nint(hproj->GetMean());
844// fPeakWidths[4]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
845 fPeakWidths[4]=(UShort_t)TMath::Nint(10.);
846
847 //other peaks
848 Int_t timepos=fPeaks[4]-2*fPeakWidths[4];
849 Int_t width=100;
850
851 for (Int_t i=3; i>=0; --i){
852 hproj->GetXaxis()->SetRangeUser(timepos-width,timepos);
853 fPeaks[i]=(UShort_t)TMath::Nint(hproj->GetMean());
854 fPeakWidths[i]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
855 width=250;
856 timepos=fPeaks[i]-width/2;
857 }
858
859 //check width and reset peak if >100
860 for (Int_t i=0; i<5; ++i){
861 if (fPeakWidths[i]>100) {
862 fPeaks[i]=0;
863 fPeakWidths[i]=0;
864 }
865 }
866
867 delete hproj;
868}
869
75d8233f 870//_____________________________________________________________________
871void AliTPCCalibCE::FindPedestal(Float_t part)
872{
ac940b58 873 //
75d8233f 874 // find pedestal and noise for the current pad. Use either database or
875 // truncated mean with part*100%
ac940b58 876 //
877 Bool_t noPedestal = kTRUE;
7fb602b1 878
879 //use pedestal database if set
ac940b58 880 if (fPedestalTPC&&fPadNoiseTPC){
75d8233f 881 //only load new pedestals if the sector has changed
ac940b58 882 if ( fCurrentSector!=fLastSector ){
883 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
884 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
885 }
75d8233f 886
ac940b58 887 if ( fPedestalROC&&fPadNoiseROC ){
888 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
889 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
890 noPedestal = kFALSE;
75d8233f 891 }
892
ac940b58 893 }
894
75d8233f 895 //if we are not running with pedestal database, or for the current sector there is no information
896 //available, calculate the pedestal and noise on the fly
ac940b58 897 if ( noPedestal ) {
898 fPadPedestal = 0;
899 fPadNoise = 0;
900 if ( fIsZeroSuppressed ) return;
901 const Int_t kPedMax = 100; //maximum pedestal value
902 Float_t max = 0;
903 Float_t maxPos = 0;
904 Int_t median = -1;
905 Int_t count0 = 0;
906 Int_t count1 = 0;
907 //
908 Float_t padSignal=0;
909 //
910 UShort_t histo[kPedMax];
911 memset(histo,0,kPedMax*sizeof(UShort_t));
bf57d87d 912
7fb602b1 913 //fill pedestal histogram
ac940b58 914 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
2963bcbf 915 padSignal = fPadSignal[i];
ac940b58 916 if (padSignal<=0) continue;
917 if (padSignal>max && i>10) {
918 max = padSignal;
919 maxPos = i;
920 }
921 if (padSignal>kPedMax-1) continue;
922 histo[int(padSignal+0.5)]++;
923 count0++;
924 }
7fb602b1 925 //find median
ac940b58 926 for (Int_t i=1; i<kPedMax; ++i){
927 if (count1<count0*0.5) median=i;
928 count1+=histo[i];
929 }
75d8233f 930 // truncated mean
ac940b58 931 //
932 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
933 //
934 for (Int_t idelta=1; idelta<10; ++idelta){
935 if (median-idelta<=0) continue;
936 if (median+idelta>kPedMax) continue;
937 if (count<part*count1){
938 count+=histo[median-idelta];
939 mean +=histo[median-idelta]*(median-idelta);
940 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
941 count+=histo[median+idelta];
942 mean +=histo[median+idelta]*(median+idelta);
943 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
944 }
945 }
946 if ( count > 0 ) {
947 mean/=count;
948 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
949 fPadPedestal = mean;
950 fPadNoise = rms;
75d8233f 951 }
ac940b58 952 }
75d8233f 953}
954//_____________________________________________________________________
ac940b58 955void AliTPCCalibCE::UpdateCETimeRef()
956{
957 // Find the time reference of the last valid CE signal in sector
958 // for irocs of the A-Side the reference of the corresponging OROC is returned
959 // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
960 if ( fLastSector == fCurrentSector ) return;
961 Int_t sector=fCurrentSector;
962 if ( sector < 18 ) sector+=36;
963 fCurrentCETimeRef=0;
964 TVectorF *vtRef = GetTMeanEvents(sector);
965 if ( !vtRef ) return;
966 Int_t vtRefSize= vtRef->GetNrows();
967 if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
968 else vtRefSize=fNevents;
969 while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
970 fCurrentCETimeRef=(*vtRef)[vtRefSize];
971 AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
972}
973//_____________________________________________________________________
75d8233f 974void AliTPCCalibCE::FindCESignal(TVectorD &param, Float_t &qSum, const TVectorF maxima)
975{
ac940b58 976 //
75d8233f 977 // Find position, signal width and height of the CE signal (last signal)
978 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
979 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
ac940b58 980 //
75d8233f 981
ac940b58 982 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
983 Int_t cemaxpos = 0;
984 Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise; // threshold for the signal sum
880c3382 985 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
986 const Int_t kCemax = fPeakIntPlus;
75d8233f 987
ac940b58 988 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
bf57d87d 989
75d8233f 990 // find maximum closest to the sector mean from the last event
ac940b58 991 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
7fb602b1 992 // get sector mean of last event
ac940b58 993 Float_t tmean = fCurrentCETimeRef;
994 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
995 minDist = tmean-maxima[imax];
996 cemaxpos = (Int_t)maxima[imax];
75d8233f 997 }
ac940b58 998 }
880c3382 999// printf("L1 phase TB: %f\n",GetL1PhaseTB());
ac940b58 1000 if (cemaxpos!=0){
2963bcbf 1001 ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
880c3382 1002 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
ac940b58 1003 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
2963bcbf 1004 Float_t signal = fPadSignal[i]-fPadPedestal;
ac940b58 1005 if (signal>0) {
1006 ceTime+=signal*(i+0.5);
1007 ceRMS +=signal*(i+0.5)*(i+0.5);
1008 ceQsum+=signal;
1009 }
1010 }
75d8233f 1011 }
ac940b58 1012 }
1013 if (ceQmax&&ceQsum>ceSumThreshold) {
1014 ceTime/=ceQsum;
1015 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
880c3382 1016 ceTime-=GetL1PhaseTB();
ac940b58 1017 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
1018 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
1019
1020 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
1021 // the pick-up signal should scale with the pad area. In addition
1022 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
1023 // ratio 2/3. The pad area we express in cm2. We normalise the signal
1024 // to the OROC signal (factor 2/3 for the IROCs).
1025 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
1026 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
1027
1028 ceQsum/=norm;
1029 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
1030 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
1031 } else {
1032 ceQmax=0;
1033 ceTime=0;
1034 ceRMS =0;
1035 ceQsum=0;
1036 }
1037 param[0] = ceQmax;
880c3382 1038 param[1] = ceTime;
ac940b58 1039 param[2] = ceRMS;
1040 qSum = ceQsum;
75d8233f 1041}
1042//_____________________________________________________________________
3cd27a08 1043Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
75d8233f 1044{
ac940b58 1045 //
2963bcbf 1046 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
1047 // and 'tplus' timebins after 'pos'
ac940b58 1048 //
1049 if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
1050 for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
1051 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
1052 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
1053 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
1054 ++iTime2;
1055 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
1056 }
1057 return kTRUE;
75d8233f 1058}
1059//_____________________________________________________________________
1060void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
1061{
ac940b58 1062 //
75d8233f 1063 // Find local maxima on the pad signal and Histogram them
ac940b58 1064 //
4c6d06dc 1065 Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal
ac940b58 1066 Int_t count = 0;
2963bcbf 1067
1068 for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
1069 if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
1070 if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
ac940b58 1071 if (count<maxima.GetNrows()){
1072 maxima.GetMatrixArray()[count++]=i;
1073 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
2963bcbf 1074 i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin fPeakDetMinus+fPeakDetPlus-1
ac940b58 1075 }
75d8233f 1076 }
ac940b58 1077 }
75d8233f 1078}
1079//_____________________________________________________________________
7fb602b1 1080void AliTPCCalibCE::ProcessPad()
75d8233f 1081{
880c3382 1082 //
1083 // Process data of current pad
1084 //
1085 FindPedestal();
1086
1087 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
7fb602b1 1088 // + central electrode and possibly post peaks from the CE signal
1089 // however if we are on a high noise pad a lot more peaks due to the noise might occur
880c3382 1090 FindLocalMaxima(maxima);
1091 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
1092
1093 UpdateCETimeRef(); // update the time refenrence for the current sector
7442bceb 1094 if ( fCurrentCETimeRef<1e-30 ) return; //return if we don't have time 0 info, eg if only one side has laser
880c3382 1095 TVectorD param(3);
1096 Float_t qSum;
1097 FindCESignal(param, qSum, maxima);
1098
1099 Double_t meanT = param[1];
1100 Double_t sigmaT = param[2];
1101
75d8233f 1102 //Fill Event T0 counter
880c3382 1103 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
1104
75d8233f 1105 //Fill Q histogram
880c3382 1106 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
1107
75d8233f 1108 //Fill RMS histogram
880c3382 1109 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
1110
1111
75d8233f 1112 //Fill debugging info
880c3382 1113 if ( GetStreamLevel()>0 ){
1114 (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
1115 (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
1116 (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
1117 }
1118
1119 ResetPad();
75d8233f 1120}
1121//_____________________________________________________________________
7fb602b1 1122void AliTPCCalibCE::EndEvent()
75d8233f 1123{
ac940b58 1124 // Process data of current pad
1125 // The Functions 'SetTimeStamp' and 'SetRunNumber' should be called
1126 // before the EndEvent function to set the event timestamp and number!!!
1127 // This is automatically done if the ProcessEvent(AliRawReader *rawReader)
1128 // function was called
78f17711 1129 if (!fProcessOld) {
1130 if (fProcessNew) ++fNevents;
1131 return;
1132 }
1133
ac940b58 1134 //check if last pad has allready been processed, if not do so
1135 if ( fMaxTimeBin>-1 ) ProcessPad();
75d8233f 1136
ac940b58 1137 AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
4c6d06dc 1138
ac940b58 1139 TVectorD param(3);
1140 TMatrixD dummy(3,3);
7fb602b1 1141// TVectorF vMeanTime(72);
1142// TVectorF vMeanQ(72);
ac940b58 1143 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
1144 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
1145
1146 //find mean time0 offset for side A and C
1147 //use only orocs due to the better statistics
1148 Double_t time0Side[2]; //time0 for side A:0 and C:1
1149 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:1
1150 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
1151 for ( Int_t iSec = 36; iSec<72; ++iSec ){
1152 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
1153 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
1154 }
1155 if ( time0SideCount[0] >0 )
1156 time0Side[0]/=time0SideCount[0];
1157 if ( time0SideCount[1] >0 )
1158 time0Side[1]/=time0SideCount[1];
75d8233f 1159 // end find time0 offset
ac940b58 1160 AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1161 Int_t nSecMeanT=0;
1162 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1163 for ( Int_t iSec = 0; iSec<72; ++iSec ){
1164 AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1165 //find median and then calculate the mean around it
1166 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
1167 if ( !hMeanT ) continue;
1168 //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1169 if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1170 hMeanT->Reset();
1171 AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1172 continue;
1173 }
1174
1175 Double_t entries = hMeanT->GetEffectiveEntries();
1176 Double_t sum = 0;
1177 Short_t *arr = hMeanT->GetArray()+1;
1178 Int_t ibin=0;
1179 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1180 sum+=arr[ibin];
1181 if ( sum>=(entries/2.) ) break;
1182 }
1183 Int_t delta = 4;
1184 Int_t firstBin = fFirstTimeBin+ibin-delta;
1185 Int_t lastBin = fFirstTimeBin+ibin+delta;
1186 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1187 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
1188 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1189
7fb602b1 1190 // check boundaries for ebye info of mean time
ac940b58 1191 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1192 Int_t vSize=vMeanTime->GetNrows();
1193 if ( vSize < fNevents+1 ){
1194 vMeanTime->ResizeTo(vSize+100);
1195 }
880c3382 1196
1197 // store mean time for the readout sides
1198 vSize=fVTime0SideA.GetNrows();
1199 if ( vSize < fNevents+1 ){
1200 fVTime0SideA.ResizeTo(vSize+100);
1201 fVTime0SideC.ResizeTo(vSize+100);
1202 }
1203 fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1204 fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
ac940b58 1205
1206 vMeanTime->GetMatrixArray()[fNevents]=median;
1207 nSecMeanT++;
1208 // end find median
1209
1210 TVectorF *vTimes = GetPadTimesEvent(iSec);
1211 if ( !vTimes ) continue; //continue if no time information for this sector is available
1212
1213 AliTPCCalROC calIrocOutliers(0);
1214 AliTPCCalROC calOrocOutliers(36);
1215
1216 // calculate mean Q of the sector
1217 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1218 vSize=vMeanQ->GetNrows();
1219 if ( vSize < fNevents+1 ){
1220 vMeanQ->ResizeTo(vSize+100);
1221 }
1222 Float_t meanQ = 0;
1223 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1224 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1225
1226 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1227 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
75d8233f 1228
1229 //set values for temporary roc calibration class
ac940b58 1230 if ( iSec < 36 ) {
1231 calIroc->SetValue(iChannel, time);
7442bceb 1232 if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
bf57d87d 1233
ac940b58 1234 } else {
1235 calOroc->SetValue(iChannel, time);
7442bceb 1236 if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
ac940b58 1237 }
75d8233f 1238
ac940b58 1239 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
c3066940 1240 // test that we really found the CE signal reliably
1241 if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1242 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
75d8233f 1243
1244
1245
1246 //------------------------------- Debug start ------------------------------
880c3382 1247 if ( GetStreamLevel()>0 ){
1248 TTreeSRedirector *streamer=GetDebugStreamer();
1249 if (streamer){
1250 Int_t row=0;
1251 Int_t pad=0;
1252 Int_t padc=0;
1253
1254 Float_t q = (*GetPadQEvent(iSec))[iChannel];
1255 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1256
1257 UInt_t channel=iChannel;
1258 Int_t sector=iSec;
1259
1260 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1261 pad = channel-fROC->GetRowIndexes(sector)[row];
1262 padc = pad-(fROC->GetNPads(sector,row)/2);
1263
75d8233f 1264// TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1265// Form("hSignalD%d.%d.%d",sector,row,pad),
1266// fLastTimeBin-fFirstTimeBin,
1267// fFirstTimeBin,fLastTimeBin);
1268// h1->SetDirectory(0);
ac940b58 1269 //
7fb602b1 1270// for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
75d8233f 1271// h1->Fill(i,fPadSignal(i));
880c3382 1272
1273 Double_t t0Sec = 0;
1274 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1275 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1276 Double_t t0Side = time0Side[(iSec/18)%2];
1277 (*streamer) << "DataPad" <<
ac940b58 1278 "Event=" << fNevents <<
1279 "RunNumber=" << fRunNumber <<
1280 "TimeStamp=" << fTimeStamp <<
1281 "Sector="<< sector <<
1282 "Row=" << row<<
1283 "Pad=" << pad <<
1284 "PadC=" << padc <<
1285 "PadSec="<< channel <<
1286 "Time0Sec=" << t0Sec <<
1287 "Time0Side=" << t0Side <<
1288 "Time=" << time <<
1289 "RMS=" << rms <<
1290 "Sum=" << q <<
1291 "MeanQ=" << meanQ <<
880c3382 1292 // "hist.=" << h1 <<
ac940b58 1293 "\n";
880c3382 1294
1295 // delete h1;
1296 }
ac940b58 1297 }
880c3382 1298 //----------------------------- Debug end ------------------------------
ac940b58 1299 }// end channel loop
1300
75d8233f 1301
2963bcbf 1302 //do fitting now only in debug mode
1303 if (GetDebugLevel()>0){
1304 TVectorD paramPol1(3);
1305 TVectorD paramPol2(6);
1306 TMatrixD matPol1(3,3);
1307 TMatrixD matPol2(6,6);
1308 Float_t chi2Pol1=0;
1309 Float_t chi2Pol2=0;
1310
1311 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1312 if ( iSec < 36 ){
1313 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1314 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1315 } else {
1316 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1317 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1318 }
1319
1320 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1321 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1322 }
1323
1324 //------------------------------- Debug start ------------------------------
1325 if ( GetStreamLevel()>0 ){
1326 TTreeSRedirector *streamer=GetDebugStreamer();
1327 if ( streamer ) {
1328 (*streamer) << "DataRoc" <<
1329// "Event=" << fEvent <<
1330 "RunNumber=" << fRunNumber <<
1331 "TimeStamp=" << fTimeStamp <<
1332 "Sector="<< iSec <<
1333 "hMeanT.=" << hMeanT <<
1334 "median=" << median <<
1335 "paramPol1.=" << &paramPol1 <<
1336 "paramPol2.=" << &paramPol2 <<
1337 "matPol1.=" << &matPol1 <<
1338 "matPol2.=" << &matPol2 <<
1339 "chi2Pol1=" << chi2Pol1 <<
1340 "chi2Pol2=" << chi2Pol2 <<
1341 "\n";
1342 }
880c3382 1343 }
ac940b58 1344 }
75d8233f 1345 //------------------------------- Debug end ------------------------------
ac940b58 1346 hMeanT->Reset();
1347 }// end sector loop
4d885988 1348 //return if no sector has a valid mean time
ac940b58 1349 if ( nSecMeanT == 0 ) return;
4d885988 1350
1351
7fb602b1 1352// fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1353// fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
ac940b58 1354 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1355 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1356 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1357 }
1358 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1359 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
75d8233f 1360
78f17711 1361 ++fNevents;
ac940b58 1362 fOldRunNumber = fRunNumber;
75d8233f 1363
ac940b58 1364 delete calIroc;
1365 delete calOroc;
1366 AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
75d8233f 1367}
1368//_____________________________________________________________________
7fb602b1 1369TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
75d8233f 1370 Int_t nbinsY, Float_t ymin, Float_t ymax,
a6e0ebfe 1371 const Char_t *type, Bool_t force)
75d8233f 1372{
1373 //
1374 // return pointer to TH2S histogram of 'type'
1375 // if force is true create a new histogram if it doesn't exist allready
1376 //
1377 if ( !force || arr->UncheckedAt(sector) )
a3b590cf 1378 return (TH2S*)arr->UncheckedAt(sector);
75d8233f 1379
7fb602b1 1380 // if we are forced and histogram doesn't exist yet create it
75d8233f 1381 // new histogram with Q calib information. One value for each pad!
a3b590cf 1382 TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
75d8233f 1383 nbinsY, ymin, ymax,
1384 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1385 hist->SetDirectory(0);
1386 arr->AddAt(hist,sector);
1387 return hist;
1388}
1389//_____________________________________________________________________
7fb602b1 1390TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
75d8233f 1391{
1392 //
1393 // return pointer to T0 histogram
1394 // if force is true create a new histogram if it doesn't exist allready
1395 //
1396 TObjArray *arr = &fHistoT0Array;
1397 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1398}
1399//_____________________________________________________________________
7fb602b1 1400TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
75d8233f 1401{
1402 //
1403 // return pointer to Q histogram
1404 // if force is true create a new histogram if it doesn't exist allready
1405 //
1406 TObjArray *arr = &fHistoQArray;
1407 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1408}
1409//_____________________________________________________________________
7fb602b1 1410TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
75d8233f 1411{
1412 //
1413 // return pointer to Q histogram
1414 // if force is true create a new histogram if it doesn't exist allready
1415 //
1416 TObjArray *arr = &fHistoRMSArray;
1417 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1418}
1419//_____________________________________________________________________
1420TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
a6e0ebfe 1421 const Char_t *type, Bool_t force)
75d8233f 1422{
1423 //
1424 // return pointer to TH1S histogram
1425 // if force is true create a new histogram if it doesn't exist allready
1426 //
1427 if ( !force || arr->UncheckedAt(sector) )
a3b590cf 1428 return (TH1S*)arr->UncheckedAt(sector);
75d8233f 1429
1430 // if we are forced and histogram doesn't yes exist create it
4d885988 1431 // new histogram with calib information. One value for each pad!
a3b590cf 1432 TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
75d8233f 1433 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1434 hist->SetDirectory(0);
1435 arr->AddAt(hist,sector);
1436 return hist;
1437}
1438//_____________________________________________________________________
1439TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1440{
1441 //
1442 // return pointer to Q histogram
1443 // if force is true create a new histogram if it doesn't exist allready
1444 //
1445 TObjArray *arr = &fHistoTmean;
1446 return GetHisto(sector, arr, "LastTmean", force);
1447}
1448//_____________________________________________________________________
3cd27a08 1449TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
75d8233f 1450{
1451 //
1452 // return pointer to Pad Info from 'arr' for the current event and sector
1453 // if force is true create it if it doesn't exist allready
1454 //
1455 if ( !force || arr->UncheckedAt(sector) )
1456 return (TVectorF*)arr->UncheckedAt(sector);
1457
7fb602b1 1458 TVectorF *vect = new TVectorF(size);
75d8233f 1459 arr->AddAt(vect,sector);
1460 return vect;
1461}
1462//_____________________________________________________________________
7fb602b1 1463TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
75d8233f 1464{
1465 //
1466 // return pointer to Pad Times Array for the current event and sector
1467 // if force is true create it if it doesn't exist allready
1468 //
1469 TObjArray *arr = &fPadTimesArrayEvent;
7fb602b1 1470 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1471}
1472//_____________________________________________________________________
7fb602b1 1473TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
75d8233f 1474{
1475 //
1476 // return pointer to Pad Q Array for the current event and sector
1477 // if force is true create it if it doesn't exist allready
1478 // for debugging purposes only
1479 //
1480
1481 TObjArray *arr = &fPadQArrayEvent;
7fb602b1 1482 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1483}
1484//_____________________________________________________________________
7fb602b1 1485TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
75d8233f 1486{
1487 //
1488 // return pointer to Pad RMS Array for the current event and sector
1489 // if force is true create it if it doesn't exist allready
1490 // for debugging purposes only
1491 //
1492 TObjArray *arr = &fPadRMSArrayEvent;
7fb602b1 1493 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
75d8233f 1494}
1495//_____________________________________________________________________
7fb602b1 1496TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
75d8233f 1497{
1498 //
1499 // return pointer to Pad RMS Array for the current event and sector
1500 // if force is true create it if it doesn't exist allready
1501 // for debugging purposes only
1502 //
1503 TObjArray *arr = &fPadPedestalArrayEvent;
7fb602b1 1504 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1505}
1506//_____________________________________________________________________
1507TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1508{
1509 //
1510 // return pointer to the EbyE info of the mean arrival time for 'sector'
1511 // if force is true create it if it doesn't exist allready
1512 //
1513 TObjArray *arr = &fTMeanArrayEvent;
1514 return GetVectSector(sector,arr,100,force);
75d8233f 1515}
1516//_____________________________________________________________________
7fb602b1 1517TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1518{
1519 //
1520 // return pointer to the EbyE info of the mean arrival time for 'sector'
1521 // if force is true create it if it doesn't exist allready
1522 //
1523 TObjArray *arr = &fQMeanArrayEvent;
1524 return GetVectSector(sector,arr,100,force);
1525}
1526//_____________________________________________________________________
3cd27a08 1527AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1528{
1529 //
1530 // return pointer to ROC Calibration
1531 // if force is true create a new histogram if it doesn't exist allready
1532 //
1533 if ( !force || arr->UncheckedAt(sector) )
1534 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1535
1536 // if we are forced and histogram doesn't yes exist create it
1537
1538 // new AliTPCCalROC for T0 information. One value for each pad!
1539 AliTPCCalROC *croc = new AliTPCCalROC(sector);
75d8233f 1540 arr->AddAt(croc,sector);
1541 return croc;
1542}
1543//_____________________________________________________________________
7fb602b1 1544AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
75d8233f 1545{
1546 //
ef7f7670 1547 // return pointer to Time 0 ROC Calibration
75d8233f 1548 // if force is true create a new histogram if it doesn't exist allready
1549 //
1550 TObjArray *arr = &fCalRocArrayT0;
1551 return GetCalRoc(sector, arr, force);
1552}
1553//_____________________________________________________________________
ef7f7670 1554AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1555{
1556 //
1557 // return pointer to the error of Time 0 ROC Calibration
1558 // if force is true create a new histogram if it doesn't exist allready
1559 //
1560 TObjArray *arr = &fCalRocArrayT0Err;
1561 return GetCalRoc(sector, arr, force);
1562}
1563//_____________________________________________________________________
7fb602b1 1564AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
75d8233f 1565{
1566 //
1567 // return pointer to T0 ROC Calibration
1568 // if force is true create a new histogram if it doesn't exist allready
1569 //
1570 TObjArray *arr = &fCalRocArrayQ;
1571 return GetCalRoc(sector, arr, force);
1572}
1573//_____________________________________________________________________
7fb602b1 1574AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
75d8233f 1575{
1576 //
1577 // return pointer to signal width ROC Calibration
1578 // if force is true create a new histogram if it doesn't exist allready
1579 //
1580 TObjArray *arr = &fCalRocArrayRMS;
1581 return GetCalRoc(sector, arr, force);
1582}
1583//_____________________________________________________________________
1584AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1585{
1586 //
1587 // return pointer to Outliers
1588 // if force is true create a new histogram if it doesn't exist allready
1589 //
1590 TObjArray *arr = &fCalRocArrayOutliers;
1591 return GetCalRoc(sector, arr, force);
1592}
1593//_____________________________________________________________________
3cd27a08 1594TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
75d8233f 1595{
1596 //
1597 // return pointer to TObjArray of fit parameters
1598 // if force is true create a new histogram if it doesn't exist allready
1599 //
1600 if ( !force || arr->UncheckedAt(sector) )
1601 return (TObjArray*)arr->UncheckedAt(sector);
1602
1603 // if we are forced and array doesn't yes exist create it
1604
1605 // new TObjArray for parameters
1606 TObjArray *newArr = new TObjArray;
1607 arr->AddAt(newArr,sector);
1608 return newArr;
1609}
1610//_____________________________________________________________________
1611TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1612{
1613 //
1614 // return pointer to TObjArray of fit parameters from plane fit
1615 // if force is true create a new histogram if it doesn't exist allready
1616 //
1617 TObjArray *arr = &fParamArrayEventPol1;
1618 return GetParamArray(sector, arr, force);
1619}
1620//_____________________________________________________________________
1621TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1622{
1623 //
1624 // return pointer to TObjArray of fit parameters from parabola fit
1625 // if force is true create a new histogram if it doesn't exist allready
1626 //
1627 TObjArray *arr = &fParamArrayEventPol2;
1628 return GetParamArray(sector, arr, force);
1629}
78f17711 1630
1631//_____________________________________________________________________
1632void AliTPCCalibCE::CreateDVhist()
1633{
1634 //
1635 // Setup the THnSparse for the drift velocity determination
1636 //
1637
1638 //HnSparse bins
1639 //roc, row, pad, timebin, timestamp
1640 TTimeStamp begin(2010,01,01,0,0,0);
1641 TTimeStamp end(2035,01,01,0,0,0);
1642 Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1643
1644 Int_t bins[kHnBinsDV] = { 72, 96, 140, 1030, nbinsTime};
1645 Double_t xmin[kHnBinsDV] = { 0., 0., 0., 0., (Double_t)begin.GetSec()};
1646 Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1647
1648 fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1649 fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1650 fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1651 fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1652 fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1653 fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
1654
1655}
1656
75d8233f 1657//_____________________________________________________________________
7fb602b1 1658void AliTPCCalibCE::ResetEvent()
75d8233f 1659{
1660 //
1661 // Reset global counters -- Should be called before each event is processed
1662 //
1663 fLastSector=-1;
1664 fCurrentSector=-1;
1665 fCurrentRow=-1;
1666 fCurrentChannel=-1;
1667
1668 ResetPad();
1669
1670 fPadTimesArrayEvent.Delete();
1671 fPadQArrayEvent.Delete();
1672 fPadRMSArrayEvent.Delete();
1673 fPadPedestalArrayEvent.Delete();
1674
7fb602b1 1675 for ( Int_t i=0; i<72; ++i ){
bf57d87d 1676 fVTime0Offset.GetMatrixArray()[i]=0;
1677 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1678 fVMeanQ.GetMatrixArray()[i]=0;
1679 fVMeanQCounter.GetMatrixArray()[i]=0;
75d8233f 1680 }
1681}
1682//_____________________________________________________________________
7fb602b1 1683void AliTPCCalibCE::ResetPad()
75d8233f 1684{
1685 //
1686 // Reset pad infos -- Should be called after a pad has been processed
1687 //
7fb602b1 1688 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
2963bcbf 1689 fPadSignal[i] = 0;
75d8233f 1690 fMaxTimeBin = -1;
1691 fMaxPadSignal = -1;
1692 fPadPedestal = -1;
1693 fPadNoise = -1;
1694}
1695//_____________________________________________________________________
7442bceb 1696void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
7fb602b1 1697{
7442bceb 1698 //
1699 // Merge ce to the current AliTPCCalibCE
1700 //
78f17711 1701 MergeBase(ce);
1702 Int_t nCEevents = ce->GetNeventsProcessed();
7442bceb 1703
78f17711 1704 if (fProcessOld&&ce->fProcessOld){
7442bceb 1705 //merge histograms
78f17711 1706 for (Int_t iSec=0; iSec<72; ++iSec){
1707 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1708 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1709 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1710
1711
1712 if ( hRefQmerge ){
1713 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1714 TH2S *hRefQ = GetHistoQ(iSec);
1715 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1716 else {
1717 TH2S *hist = new TH2S(*hRefQmerge);
1718 hist->SetDirectory(0);
1719 fHistoQArray.AddAt(hist, iSec);
1720 }
1721 hRefQmerge->SetDirectory(dir);
7442bceb 1722 }
78f17711 1723 if ( hRefT0merge ){
1724 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1725 TH2S *hRefT0 = GetHistoT0(iSec);
1726 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1727 else {
1728 TH2S *hist = new TH2S(*hRefT0merge);
1729 hist->SetDirectory(0);
1730 fHistoT0Array.AddAt(hist, iSec);
1731 }
1732 hRefT0merge->SetDirectory(dir);
7442bceb 1733 }
78f17711 1734 if ( hRefRMSmerge ){
1735 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1736 TH2S *hRefRMS = GetHistoRMS(iSec);
1737 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1738 else {
1739 TH2S *hist = new TH2S(*hRefRMSmerge);
1740 hist->SetDirectory(0);
1741 fHistoRMSArray.AddAt(hist, iSec);
1742 }
1743 hRefRMSmerge->SetDirectory(dir);
7442bceb 1744 }
78f17711 1745
7442bceb 1746 }
1747
7fb602b1 1748 // merge time information
7442bceb 1749
7442bceb 1750
78f17711 1751 for (Int_t iSec=0; iSec<72; ++iSec){
1752 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1753 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1754 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1755 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1756
1757 TObjArray *arrPol1 = 0x0;
1758 TObjArray *arrPol2 = 0x0;
1759 TVectorF *vMeanTime = 0x0;
1760 TVectorF *vMeanQ = 0x0;
1761
7442bceb 1762 //resize arrays
78f17711 1763 if ( arrPol1CE && arrPol2CE ){
1764 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1765 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1766 arrPol1->Expand(fNevents+nCEevents);
1767 arrPol2->Expand(fNevents+nCEevents);
1768 }
1769 if ( vMeanTimeCE && vMeanQCE ){
1770 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1771 vMeanQ = GetQMeanEvents(iSec,kTRUE);
1772 vMeanTime->ResizeTo(fNevents+nCEevents);
1773 vMeanQ->ResizeTo(fNevents+nCEevents);
1774 }
1775
1776 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1777 if ( arrPol1CE && arrPol2CE ){
1778 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1779 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1780 if ( paramPol1 && paramPol2 ){
1781 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1782 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1783 }
1784 }
1785 if ( vMeanTimeCE && vMeanQCE ){
1786 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1787 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1788 }
1789 }
7442bceb 1790 }
1791
1792
78f17711 1793
1794 const TVectorD& eventTimes = ce->fVEventTime;
1795 const TVectorD& eventIds = ce->fVEventNumber;
1796 const TVectorF& time0SideA = ce->fVTime0SideA;
1797 const TVectorF& time0SideC = ce->fVTime0SideC;
1798 fVEventTime.ResizeTo(fNevents+nCEevents);
1799 fVEventNumber.ResizeTo(fNevents+nCEevents);
1800 fVTime0SideA.ResizeTo(fNevents+nCEevents);
1801 fVTime0SideC.ResizeTo(fNevents+nCEevents);
1802
7442bceb 1803 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
78f17711 1804 Double_t evTime = eventTimes.GetMatrixArray()[iEvent];
1805 Double_t evId = eventIds.GetMatrixArray()[iEvent];
1806 Float_t t0SideA = time0SideA.GetMatrixArray()[iEvent];
1807 Float_t t0SideC = time0SideC.GetMatrixArray()[iEvent];
1808
1809 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1810 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1811 fVTime0SideA.GetMatrixArray()[fNevents+iEvent] = t0SideA;
1812 fVTime0SideC.GetMatrixArray()[fNevents+iEvent] = t0SideC;
1813 }
1814 }
1815
1816 if (fProcessNew&&ce->fProcessNew) {
1817 if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1818 AliError("Number of bursts in the instances to merge are different. No merging done!");
1819 } else {
1820 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1821 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1822 THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1823 if (h && hce) h->Add(hce);
1824 else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
7442bceb 1825 }
78f17711 1826 //TODO: What to do with fTimeBursts???
7fb602b1 1827 }
7442bceb 1828 }
1829
7442bceb 1830 fNevents+=nCEevents; //increase event counter
1831}
7fb602b1 1832
7442bceb 1833//_____________________________________________________________________
1834Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1835{
1836 //
1837 // Merge all objects of this type in list
1838 //
7fb602b1 1839
7442bceb 1840 Long64_t nmerged=1;
7fb602b1 1841
7442bceb 1842 TIter next(list);
1843 AliTPCCalibCE *ce=0;
1844 TObject *o=0;
7fb602b1 1845
7442bceb 1846 while ( (o=next()) ){
1847 ce=dynamic_cast<AliTPCCalibCE*>(o);
1848 if (ce){
1849 Merge(ce);
1850 ++nmerged;
7fb602b1 1851 }
7442bceb 1852 }
7fb602b1 1853
7442bceb 1854 return nmerged;
7fb602b1 1855}
7442bceb 1856
7fb602b1 1857//_____________________________________________________________________
1858TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
75d8233f 1859{
ac940b58 1860 //
1861 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
880c3382 1862 // or side (-1: A-Side, -2: C-Side)
ac940b58 1863 // xVariable: 0-event time, 1-event id, 2-internal event counter
1864 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1865 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1866 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1867 // not used for mean time and mean Q )
1868 // for an example see class description at the beginning
1869 //
75d8233f 1870
ac940b58 1871 TVectorD *xVar = 0x0;
1872 TObjArray *aType = 0x0;
1873 Int_t npoints=0;
75d8233f 1874
96bf9029 1875 // sanity checks
1876 if ( (sector<-2) || (sector>71) ) return 0x0; //sector outside valid range
1877 if ( (xVariable<0) || (xVariable>2) ) return 0x0; //invalid x-variable
1878 if ( (fitType<0) || (fitType>3) ) return 0x0; //invalid fit type
1879 if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1880 if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1881 if ( sector<0 && fitType!=2) return 0x0; //for side wise information only mean time is available
880c3382 1882
1883 if (sector>=0){
1884 if ( fitType==0 ){
1885 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1886 aType = &fParamArrayEventPol1;
1887 if ( aType->At(sector)==0x0 ) return 0x0;
1888 }
1889 else if ( fitType==1 ){
1890 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1891 aType = &fParamArrayEventPol2;
1892 if ( aType->At(sector)==0x0 ) return 0x0;
1893 }
75d8233f 1894
880c3382 1895 }
ac940b58 1896 if ( xVariable == 0 ) xVar = &fVEventTime;
1897 if ( xVariable == 1 ) xVar = &fVEventNumber;
1898 if ( xVariable == 2 ) {
1899 xVar = new TVectorD(fNevents);
1900 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1901 }
880c3382 1902
a3b590cf 1903 Double_t *x = new Double_t[fNevents];
1904 Double_t *y = new Double_t[fNevents];
1905
ac940b58 1906 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1907 if ( fitType<2 ){
1908 TObjArray *events = (TObjArray*)(aType->At(sector));
1909 if ( events->GetSize()<=ievent ) break;
1910 TVectorD *v = (TVectorD*)(events->At(ievent));
1911 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1912 } else if (fitType == 2) {
1913 Double_t xValue=(*xVar)[ievent];
880c3382 1914 Double_t yValue=0;
1915 if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1916 else if (sector==-1) yValue=fVTime0SideA(ievent);
1917 else if (sector==-2) yValue=fVTime0SideC(ievent);
ac940b58 1918 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1919 }else if (fitType == 3) {
1920 Double_t xValue=(*xVar)[ievent];
1921 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1922 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
75d8233f 1923 }
ac940b58 1924 }
75d8233f 1925
ac940b58 1926 TGraph *gr = new TGraph(npoints);
75d8233f 1927 //sort xVariable increasing
ac940b58 1928 Int_t *sortIndex = new Int_t[npoints];
a3b590cf 1929 TMath::Sort(npoints,x,sortIndex, kFALSE);
ac940b58 1930 for (Int_t i=0;i<npoints;++i){
1931 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1932 }
75d8233f 1933
1934
ac940b58 1935 if ( xVariable == 2 ) delete xVar;
4ce766eb 1936 delete [] x;
1937 delete [] y;
1938 delete [] sortIndex;
ac940b58 1939 return gr;
75d8233f 1940}
1941//_____________________________________________________________________
1942void AliTPCCalibCE::Analyse()
1943{
880c3382 1944 //
1945 // Calculate calibration constants
1946 //
78f17711 1947
1948 if (fProcessOld){
1949 TVectorD paramQ(3);
1950 TVectorD paramT0(3);
1951 TVectorD paramRMS(3);
1952 TMatrixD dummy(3,3);
1953
1954 Float_t channelCounter=0;
1955 fMeanT0rms=0;
1956 fMeanQrms=0;
1957 fMeanRMSrms=0;
1958
1959 for (Int_t iSec=0; iSec<72; ++iSec){
1960 TH2S *hT0 = GetHistoT0(iSec);
1961 if (!hT0 ) continue;
1962
1963 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1964 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1965 AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1966 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1967 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1968
1969 TH2S *hQ = GetHistoQ(iSec);
1970 TH2S *hRMS = GetHistoRMS(iSec);
1971
1972 Short_t *arrayhQ = hQ->GetArray();
1973 Short_t *arrayhT0 = hT0->GetArray();
1974 Short_t *arrayhRMS = hRMS->GetArray();
1975
1976 UInt_t nChannels = fROC->GetNChannels(iSec);
1977
1978 //debug
1979 Int_t row=0;
1980 Int_t pad=0;
1981 Int_t padc=0;
1982 //! debug
1983
1984 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1985
1986
1987 Float_t cogTime0 = -1000;
1988 Float_t cogQ = -1000;
1989 Float_t cogRMS = -1000;
1990 Float_t cogOut = 0;
1991 Float_t rms = 0;
1992 Float_t rmsT0 = 0;
1993
1994
1995 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1996 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1997 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
1998
1999 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2000 fMeanQrms+=rms;
2001 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2002 fMeanT0rms+=rmsT0;
2003 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2004 fMeanRMSrms+=rms;
2005 channelCounter++;
2006
2007 /*
2008 //outlier specifications
2009 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2010 cogOut = 1;
2011 cogTime0 = 0;
2012 cogQ = 0;
2013 cogRMS = 0;
2014 }
2015 */
2016 rocQ->SetValue(iChannel, cogQ*cogQ);
2017 rocT0->SetValue(iChannel, cogTime0);
2018 rocT0Err->SetValue(iChannel, rmsT0);
2019 rocRMS->SetValue(iChannel, cogRMS);
2020 rocOut->SetValue(iChannel, cogOut);
2021
2022
2023 //debug
2024 if ( GetStreamLevel() > 0 ){
2025 TTreeSRedirector *streamer=GetDebugStreamer();
2026 if ( streamer ) {
2027
2028 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2029 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2030 padc = pad-(fROC->GetNPads(iSec,row)/2);
2031
2032 (*streamer) << "DataEnd" <<
2033 "Sector=" << iSec <<
2034 "Pad=" << pad <<
2035 "PadC=" << padc <<
2036 "Row=" << row <<
2037 "PadSec=" << iChannel <<
2038 "Q=" << cogQ <<
2039 "T0=" << cogTime0 <<
2040 "RMS=" << cogRMS <<
2041 "\n";
2042 }
2043 }
2044 //! debug
2045
2046 }
2047
2048 }
2049 if ( channelCounter>0 ){
2050 fMeanT0rms/=channelCounter;
2051 fMeanQrms/=channelCounter;
2052 fMeanRMSrms/=channelCounter;
2053 }
2054 // if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2055 // delete fDebugStreamer;
2056 // fDebugStreamer = 0x0;
2057 fVEventTime.ResizeTo(fNevents);
2058 fVEventNumber.ResizeTo(fNevents);
2059 fVTime0SideA.ResizeTo(fNevents);
2060 fVTime0SideC.ResizeTo(fNevents);
2061 }
2062
2063 if (fProcessNew&&fAnalyseNew){
2064 AnalyseTrack();
2065 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2066 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2067 h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2068 }
2069 }
2070}
2071
2072
2073
2074
2075//
2076// New functions that also use the laser tracks
2077//
2078
2079
2080
2081//_____________________________________________________________________
2082void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2083{
2084 //
2085 //Find the local maximums for the projections to each axis
2086 //
880c3382 2087
78f17711 2088 //find laser layer positoins
2089 fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2090 FindLaserLayers();
2091 THnSparse *hProj=fHnDrift;
2092 Double_t posCE[4]={0.,0.,0.,0.};
2093 Double_t widthCE[4]={0.,0.,0.,0.};
880c3382 2094
78f17711 2095// if(fPeaks[4]!=0){
2096 // find central electrode position once more, separately for IROC, OROC, A-, C-Side
880c3382 2097
78f17711 2098 for (Int_t i=0; i<4; ++i){
2099 hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2100 TH1 *h=fHnDrift->Projection(3);
2101 h->GetXaxis()->SetRangeUser(fPeaks[4]-fPeakWidths[4],fPeaks[4]+fPeakWidths[4]);
2102 Int_t nbinMax=h->GetMaximumBin();
2103 Double_t maximum=h->GetMaximum();
2104// Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2105// if (nbinMax<700||maximum<maxExpected) continue;
2106 Double_t xbinMax=h->GetBinCenter(nbinMax);
2107 TF1 fgaus("gaus","gaus",xbinMax-10,xbinMax+10);
2108 fgaus.SetParameters(maximum,xbinMax,2);
2109 fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2110 fgaus.SetParLimits(2,0.2,4.);
2111 h->Fit(&fgaus,"RQN");
2112// Double_t deltaX=4*fgaus.GetParameter(2);
2113// xbinMax=fgaus.GetParameter(1);
2114 delete h;
2115 posCE[i]=fgaus.GetParameter(1);
2116 widthCE[i]=4*fgaus.GetParameter(2);
2117 hProj->GetAxis(0)->SetRangeUser(0,72);
2118 }
2119// }
2120 //Current drift velocity
2121 Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2122// cout<<"vdrift="<<vdrift<<endl;
2123
2124 AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2125 //loop over all entries in the histogram
2126 Int_t coord[5];
2127 for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2128 //get entry position and content
2129 Double_t adc=hProj->GetBinContent(ichunk,coord);
2130
880c3382 2131
78f17711 2132 Int_t sector = coord[0]-1;
2133 Int_t row = coord[1]-1;
2134 Int_t pad = coord[2]-1;
2135 Int_t timeBin= coord[3]-1;
2136 Double_t time = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2137 Int_t side = (sector/18)%2;
2138// return;
2139// fPeaks[4]=(UInt_t)posCE[sector/18];
2140// fPeakWidths[4]=(UInt_t)widthCE[sector/18];
880c3382 2141
78f17711 2142 //cuts
2143 if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2144 if (adc < 5 ) continue;
2145 if (IsEdgePad(sector,row,pad)) continue;
2146// if (!IsPeakInRange(timeBin)) continue;
2147// if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2148// if (isCE&&(sector!=0)) continue;
880c3382 2149
78f17711 2150 Int_t padmin=-2, padmax=2;
2151 Int_t timemin=-2, timemax=2;
2152 Int_t minsumperpad=2;
2153 //CE or laser tracks
2154 Bool_t isCE=kFALSE;
2155 if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2156 isCE=kTRUE;
2157 padmin=0;
2158 padmax=0;
2159 timemin=-3;
2160 timemax=7;
2161 }
880c3382 2162
78f17711 2163 //
2164 // Find local maximum and cogs
2165 //
2166 Bool_t isMaximum=kTRUE;
2167 Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2168 Double_t cogY=0, rmsY=0;
2169 Int_t npart=0;
880c3382 2170
78f17711 2171 // for position calculation use
2172 for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2173 Float_t lxyz[3];
2174 fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2175
2176 for(Int_t itime=timemin;itime<=timemax;++itime){
2177
2178 Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2179 Double_t val=hProj->GetBinContent(a);
2180 totalmass+=val;
2181
2182 tbcm +=(timeBin+itime)*val;
2183 padcm+=(pad+ipad)*val;
2184 cogY +=lxyz[1]*val;
2185
2186 rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2187 rmspad+=(pad+ipad)*(pad+ipad)*val;
2188 rmsY +=lxyz[1]*lxyz[1]*val;
2189
2190 if (val>0) ++npart;
2191 if (val>adc) {
2192 isMaximum=kFALSE;
2193 break;
2194 }
2195 }
2196 if (!isMaximum) break;
2197 }
2198
2199 if (!isMaximum||!npart) continue;
2200 if (totalmass<npart*minsumperpad) continue;
2201 if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2202 //for CE we want only one pad by construction
880c3382 2203
78f17711 2204 tbcm/=totalmass;
2205 padcm/=totalmass;
2206 cogY/=totalmass;
2207
2208 rmstb/=totalmass;
2209 rmspad/=totalmass;
2210 rmsY/=totalmass;
2211
2212 rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2213 rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2214 rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2215
2216 Int_t cog=TMath::Nint(padcm);
2217
2218 // timebin --> z position
2219 Float_t zlength=fParam->GetZLength(side);
2220// Float_t timePos=tbcm+(1000-fPeaks[4])
2221 // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2222 Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2223
2224 // local to global transformation--> x and y positions
2225 Float_t padlxyz[3];
2226 fROC->GetPositionLocal(sector,row,pad,padlxyz);
2227
2228 Double_t gxyz[3]={padlxyz[0],cogY,gz};
2229 Double_t lxyz[3]={padlxyz[0],cogY,gz};
2230 Double_t igxyz[3]={0,0,0};
2231 AliTPCTransform t1;
2232 t1.RotatedGlobal2Global(sector,gxyz);
2233
2234 Double_t mindist=0;
2235 Int_t trackID=-1;
2236 Int_t trackID2=-1;
2237
2238 //find track id and index of the position in the track (row)
2239 Int_t index=0;
2240 if (!isCE){
2241 index=row+(sector>35)*fROC->GetNRows(0);
2242 trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2243 } else {
2244 trackID=336+((sector/18)%2);
2245 index= fROC->GetRowIndexes(sector)[row]+pad; // global pad position in sector
2246 if (sector<36) {
2247 index+=(sector%18)*fROC->GetNChannels(sector);
2248 } else {
2249 index+=18*fROC->GetNChannels(0);
2250 index+=(sector%18)*fROC->GetNChannels(sector);
2251 }
2252 //TODO: find out about the multiple peaks in the CE
2253// mindist=TMath::Abs(fPeaks[4]-tbcm);
2254 mindist=1.;
2255 }
2256
2257 // fill track vectors
2258 if (trackID>0){
2259 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2260 Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2261
880c3382 2262
78f17711 2263// travel time effect of light includes
880c3382 2264
78f17711 2265 Double_t raylength=ltr->GetRayLength();
2266 Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2267 ltr->GetXYZ(globmir);
2268 if(trackID<336){
2269 if(side==0){
2270 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2271 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
a7307087 2272 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
78f17711 2273 }
2274 else {
2275 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2276 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
a7307087 2277 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
78f17711 2278 }
2279 }
880c3382 2280
78f17711 2281 if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2282 ltr->fVecSec->GetMatrixArray()[index]=sector;
2283 ltr->fVecP2->GetMatrixArray()[index]=0;
2284 ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2285 ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2286 ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2287 ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2288 ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2289 ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2290 ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2291// ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2292 }
2293 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2294 ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2295 igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2296 igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2297 igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2298 }
2299
2300
2301 if (fStreamLevel>4){
2302 (*GetDebugStreamer()) << "clusters" <<
2303 "run=" << fRunNumber <<
2304 "timestamp=" << timestamp <<
2305 "burst=" << burst <<
2306 "side=" << side <<
2307 "sec=" << sector <<
2308 "row=" << row <<
2309 "pad=" << pad <<
2310 "padCog=" << cog <<
2311 "timebin=" << timeBin <<
2312 "cogCE=" << posCE[sector/18] <<
2313 "withCE=" << widthCE[sector/18] <<
2314 "index=" << index <<
2315
2316 "padcm=" << padcm <<
2317 "rmspad=" << rmspad <<
2318
2319 "cogtb=" << tbcm <<
2320 "rmstb=" << rmstb <<
2321
2322 "npad=" << npart <<
2323
2324 "lx=" << padlxyz[0]<<
2325 "ly=" << cogY <<
2326 "lypad=" << padlxyz[1]<<
2327 "rmsY=" << rmsY <<
2328
2329 "gx=" << gxyz[0] <<
2330 "gy=" << gxyz[1] <<
2331 "gz=" << gxyz[2] <<
2332
2333 "igx=" << igxyz[0] <<
2334 "igy=" << igxyz[1] <<
2335 "igz=" << igxyz[2] <<
2336
2337 "mind=" << mindist <<
2338 "max=" << adc <<
2339 "trackid=" << trackID <<
2340 "trackid2=" << trackID2 <<
2341 "npart=" << npart <<
2342 "\n";
2343 } // end stream levelmgz.fElements
2344
2345 }
2346
2347}
2348
2349//_____________________________________________________________________
2350void AliTPCCalibCE::AnalyseTrack()
2351{
2352 //
2353 // Analyse the tracks
2354 //
2355
2356
2357 AliTPCLaserTrack::LoadTracks();
2358// AliTPCParam *param=0x0;
2359// //cdb run number
2360// AliCDBManager *man=AliCDBManager::Instance();
2361// if (man->GetDefaultStorage()){
2362// AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2363// if (entry){
2364// entry->SetOwner(kTRUE);
2365// param = (AliTPCParam*)(entry->GetObject()->Clone());
2366// }
2367// }
2368// if (param){
2369// if (fParam) delete fParam;
2370// fParam=param;
2371// } else {
2372// AliError("Could not get updated AliTPCParam from OCDB!!!");
2373// }
2374
2375 //Measured and ideal laser tracks
2376 TObjArray* arrMeasured = SetupMeasured();
2377 TObjArray* arrIdeal = AliTPCLaserTrack::GetTracks();
2378 AddCEtoIdeal(arrIdeal);
2379
2380 //find bursts and loop over them
2381 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2382 Double_t timestamp=fTimeBursts[iburst];
2383 AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2384 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2385 if (!fHnDrift) continue;
2386 UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2387 if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2388 fBinsLastAna[iburst]=entries;
2389
2390 for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2391// if (iburst==0) FindLaserLayers();
2392
2393 //reset laser tracks
2394 ResetMeasured(arrMeasured);
2395
2396 //find clusters and associate to the tracks
2397 FindLocalMaxima(arrMeasured, timestamp, iburst);
2398
2399 //calculate drift velocity
2400 CalculateDV(arrIdeal,arrMeasured,iburst);
2401
2402 //Dump information to file if requested
2403 if (fStreamLevel>2){
2404 printf("make tree\n");
2405 //laser track information
880c3382 2406
78f17711 2407 for (Int_t itrack=0; itrack<338; ++itrack){
2408 TObject *iltr=arrIdeal->UncheckedAt(itrack);
2409 TObject *mltr=arrMeasured->UncheckedAt(itrack);
2410 (*GetDebugStreamer()) << "tracks" <<
2411 "run=" << fRunNumber <<
2412 "time=" << timestamp <<
2413 "burst="<< iburst <<
2414 "iltr.=" << iltr <<
2415 "mltr.=" << mltr <<
2416 "\n";
2417 }
2418 }
2419 }
2420 if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2421}
2422
2423//_____________________________________________________________________
2424Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2425 const Double_t *peakposloc, Int_t &itrackMin2)
2426{
2427 //
2428 // Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2429 //
2430
2431
2432 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2433 TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2434 TVector3 vDir;
2435 TVector3 vSt;
2436
2437 Int_t firstbeam=0;
2438 Int_t lastbeam=336/2;
2439 if ( (sector/18)%2 ) {
2440 firstbeam=336/2;
2441 lastbeam=336;
2442 }
2443
2444 mindist=1000000;
2445 Int_t itrackMin=-1;
2446 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2447 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2448// if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2449 vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2450 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2451 if (TMath::Abs(deltaZ)>40) continue;
2452 vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2453 vSt.RotateZ(ltr->GetAlpha());
2454 vDir.RotateZ(ltr->GetAlpha());
2455
2456 Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2457
2458 if (dist<mindist){
2459 mindist=dist;
2460 itrackMin=itrack;
2461 }
2462
2463 }
2464 itrackMin2=-1;
2465 Float_t mindist2=10;
2466 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2467 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2468 if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2469
2470 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2471 if (TMath::Abs(deltaZ)>40) continue;
2472
2473 Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2474 if (dist>1) continue;
2475
2476 if (dist<mindist2){
2477 mindist2=dist;
2478 itrackMin2=itrack;
2479 }
2480 }
2481 mindist=mindist2;
2482 return itrackMin2;
2483
2484}
2485
2486//_____________________________________________________________________
2487Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2488{
2489 //
2490 // return true if pad is on the edge of a row
2491 //
2492 Int_t edge1 = 0;
2493 if ( pad == edge1 ) return kTRUE;
2494 Int_t edge2 = fROC->GetNPads(sector,row)-1;
2495 if ( pad == edge2 ) return kTRUE;
2496
2497 return kFALSE;
2498}
2499
2500//_____________________________________________________________________
2501TObjArray* AliTPCCalibCE::SetupMeasured()
2502{
2503 //
2504 // setup array of measured laser tracks and CE
2505 //
2506
2507 TObjArray *arrIdeal = AliTPCLaserTrack::GetTracks();
2508 TObjArray *arrMeasured = new TObjArray(338);
2509 arrMeasured->SetOwner();
2510 for(Int_t itrack=0;itrack<336;itrack++){
2511 arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2512 }
2513
2514 //"tracks" for CE
2515 AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2516 ltrce->SetId(336);
2517 ltrce->SetSide(0);
2518 ltrce->fVecSec=new TVectorD(557568/2);
2519 ltrce->fVecP2=new TVectorD(557568/2);
2520 ltrce->fVecPhi=new TVectorD(557568/2);
2521 ltrce->fVecGX=new TVectorD(557568/2);
2522 ltrce->fVecGY=new TVectorD(557568/2);
2523 ltrce->fVecGZ=new TVectorD(557568/2);
2524 ltrce->fVecLX=new TVectorD(557568/2);
2525 ltrce->fVecLY=new TVectorD(557568/2);
2526 ltrce->fVecLZ=new TVectorD(557568/2);
2527
2528 arrMeasured->AddAt(ltrce,336); //CE A-Side
2529
2530 ltrce=new AliTPCLaserTrack;
2531 ltrce->SetId(337);
2532 ltrce->SetSide(1);
2533 ltrce->fVecSec=new TVectorD(557568/2);
2534 ltrce->fVecP2=new TVectorD(557568/2);
2535 ltrce->fVecPhi=new TVectorD(557568/2);
2536 ltrce->fVecGX=new TVectorD(557568/2);
2537 ltrce->fVecGY=new TVectorD(557568/2);
2538 ltrce->fVecGZ=new TVectorD(557568/2);
2539 ltrce->fVecLX=new TVectorD(557568/2);
2540 ltrce->fVecLY=new TVectorD(557568/2);
2541 ltrce->fVecLZ=new TVectorD(557568/2);
2542 arrMeasured->AddAt(ltrce,337); //CE C-Side
2543
2544 return arrMeasured;
2545}
2546
2547//_____________________________________________________________________
2548void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2549{
2550 //
2551 // reset array of measured laser tracks and CE
2552 //
2553 for(Int_t itrack=0;itrack<338;itrack++){
2554 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2555 ltr->fVecSec->Zero();
2556 ltr->fVecP2->Zero();
2557 ltr->fVecPhi->Zero();
2558 ltr->fVecGX->Zero();
2559 ltr->fVecGY->Zero();
2560 ltr->fVecGZ->Zero();
2561 ltr->fVecLX->Zero();
2562 ltr->fVecLY->Zero();
2563 ltr->fVecLZ->Zero();
2564 }
2565}
2566
2567//_____________________________________________________________________
2568void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2569{
2570 //
2571 // Add ideal CE positions to the ideal track data
2572 //
2573
2574 arr->Expand(338);
2575 //"tracks" for CE
2576 AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2577 ltrceA->SetId(336);
2578 ltrceA->SetSide(0);
2579 ltrceA->fVecSec=new TVectorD(557568/2);
2580 ltrceA->fVecP2=new TVectorD(557568/2);
2581 ltrceA->fVecPhi=new TVectorD(557568/2);
2582 ltrceA->fVecGX=new TVectorD(557568/2);
2583 ltrceA->fVecGY=new TVectorD(557568/2);
2584 ltrceA->fVecGZ=new TVectorD(557568/2);
2585 ltrceA->fVecLX=new TVectorD(557568/2);
2586 ltrceA->fVecLY=new TVectorD(557568/2);
2587 ltrceA->fVecLZ=new TVectorD(557568/2);
2588 arr->AddAt(ltrceA,336); //CE A-Side
2589
2590 AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2591 ltrceC->SetId(337);
2592 ltrceC->SetSide(1);
2593 ltrceC->fVecSec=new TVectorD(557568/2);
2594 ltrceC->fVecP2=new TVectorD(557568/2);
2595 ltrceC->fVecPhi=new TVectorD(557568/2);
2596 ltrceC->fVecGX=new TVectorD(557568/2);
2597 ltrceC->fVecGY=new TVectorD(557568/2);
2598 ltrceC->fVecGZ=new TVectorD(557568/2);
2599 ltrceC->fVecLX=new TVectorD(557568/2);
2600 ltrceC->fVecLY=new TVectorD(557568/2);
2601 ltrceC->fVecLZ=new TVectorD(557568/2);
2602 arr->AddAt(ltrceC,337); //CE C-Side
2603
2604 //Calculate ideal positoins
2605 Float_t gxyz[3];
2606 Float_t lxyz[3];
2607 Int_t channelSideA=0;
2608 Int_t channelSideC=0;
2609 Int_t channelSide=0;
2610 AliTPCLaserTrack *ltrce=0x0;
2611
2612 for (Int_t isector=0; isector<72; ++isector){
2613 Int_t side=((isector/18)%2);
2614 for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2615 for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2616 fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2617 fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2618 if (side==0) {
2619 ltrce=ltrceA;
2620 channelSide=channelSideA;
2621 } else {
2622 ltrce=ltrceC;
2623 channelSide=channelSideC;
2624 }
2625
2626 ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2627 ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2628 ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2629 ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2630 ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2631// ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2632 ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2633 ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2634// ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2635
2636 if (side==0){
2637 ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2638 ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2639 ++channelSideA;
2640 }
2641 else {
2642 ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2643 ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2644 ++channelSideC;
2645 }
2646 }
2647 }
2648 }
2649
2650
2651}
2652
2653//_____________________________________________________________________
2654void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2655{
2656 //
2657 // calculate the drift velocity from the reconstructed clusters associated
2658 // to the ideal laser tracks
2659 // use two different fit scenarios: Separate fit for A- and C-Side
2660 // Common fit for A- and C-Side
2661 //
2662
2663 if (!fArrFitGraphs){
2664 fArrFitGraphs=new TObjArray;
2665 }
2666
2667// static TLinearFitter fdriftA(5,"hyp4");
2668// static TLinearFitter fdriftC(5,"hyp4");
2669// static TLinearFitter fdriftAC(6,"hyp5");
2670 Double_t timestamp=fTimeBursts[burst];
2671
2672 static TLinearFitter fdriftA(4,"hyp3");
2673 static TLinearFitter fdriftC(4,"hyp3");
2674 static TLinearFitter fdriftAC(5,"hyp4");
2675 TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2676
2677 Float_t chi2A = 10;
2678 Float_t chi2C = 10;
2679 Float_t chi2AC = 10;
2680 Int_t npointsA=0;
2681 Int_t npointsC=0;
2682 Int_t npointsAC=0;
2683
2684 Double_t minres[3]={20.,1,0.8};
2685 //----
2686 for(Int_t i=0;i<3;i++){
2687
2688 fdriftA.ClearPoints();
2689 fdriftC.ClearPoints();
2690 fdriftAC.ClearPoints();
2691
2692 chi2A = 10;
2693 chi2C = 10;
2694 chi2AC = 10;
2695 npointsA=0;
2696 npointsC=0;
2697 npointsAC=0;
2698
2699 for (Int_t itrack=0; itrack<338; ++itrack){
2700 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2701 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2702
2703 //-- Exclude the tracks which has the biggest inclanation angle
2704 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2705 Int_t clustercounter=0;
2706 Int_t indexMax=159;
880c3382 2707
78f17711 2708 //-- exclude the low intensity tracks
880c3382 2709
78f17711 2710 for (Int_t index=0; index<indexMax; ++index){
2711
2712 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2713 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2714 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2715
2716 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
880c3382 2717 }
78f17711 2718 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2719 clustercounter=0;
2720
880c3382 2721
78f17711 2722 //-- drift length
2723 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
880c3382 2724
78f17711 2725 if (itrack>335) indexMax=557568/2;
2726 for (Int_t index=0; index<indexMax; ++index){
2727 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2728 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2729 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2730 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
880c3382 2731
78f17711 2732 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2733 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2734 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2735 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
880c3382 2736
78f17711 2737 //cut if no track info available
2738 if (iltr->GetBundle()==0) continue;
2739 if (iR<133||mR<133) continue;
2740 if(mltr->fVecP2->GetMatrixArray()[index]>minres[i]) continue;
2741
2742 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2743 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2744
2745 //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2746 Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2747
2748 if (iltr->GetSide()==0){
2749 fdriftA.AddPoint(xxx,mdrift,1);
2750 }else{
2751 fdriftC.AddPoint(xxx,mdrift,1);
880c3382 2752 }
78f17711 2753// Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2754 Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2755 fdriftAC.AddPoint(xxx2,mdrift,1);
2756
2757 }//end index loop
2758 }//end laser track loop
2759
2760 //perform fit
2761 fdriftA.Eval();
2762 fdriftC.Eval();
2763 fdriftAC.Eval();
2764
2765
2766
2767 //get fit values
2768 fdriftA.GetParameters(fitA);
2769 fdriftC.GetParameters(fitC);
2770 fdriftAC.GetParameters(fitAC);
2771
2772 //Parameters: 0 linear offset
2773 // 1 mean drift velocity correction factor
2774 // 2 relative global y gradient
2775 // 3 radial deformation
2776 // 4 IROC/OROC offset
2777
2778// FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2779
2780 for (Int_t itrack=0; itrack<338; ++itrack){
2781 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2782 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2783
2784 //-- Exclude the tracks which has the biggest inclanation angle
2785 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2786 Int_t clustercounter=0;
2787 Int_t indexMax=159;
2788
2789 //-- exclude the low intensity tracks
2790
2791 for (Int_t index=0; index<indexMax; ++index){
2792 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2793 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2794 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2795 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
880c3382 2796 }
78f17711 2797 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2798 clustercounter=0;
2799
2800 //-- drift length
2801 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2802
2803 if (itrack>335) indexMax=557568/2;
2804 for (Int_t index=0; index<indexMax; ++index){
2805 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2806 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2807 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2808 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2809
2810 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2811 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2812 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2813 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2814
2815 //cut if no track info available
2816 if (iR<60||mR<60) continue;
2817
2818 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2819 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2820
2821 TVectorD *v=&fitA;
2822 if (iltr->GetSide()==1) v=&fitC;
2823// Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR)+(*v)[4]*( iltr->fVecSec->GetMatrixArray()[index]>35);
2824 Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2825
2826 mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2827
2828 }
2829 }
2830
2831 fitA.ResizeTo(7);
2832 fitC.ResizeTo(7);
2833 fitAC.ResizeTo(8);
2834
2835//set statistics values
2836
2837 npointsA= fdriftA.GetNpoints();
2838 if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2839 fitA[5]=npointsA;
2840 fitA[6]=chi2A;
2841
2842 npointsC= fdriftC.GetNpoints();
2843 if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2844 fitC[5]=npointsC;
2845 fitC[6]=chi2C;
2846
2847 npointsAC= fdriftAC.GetNpoints();
2848 if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2849 fitAC[5]=npointsAC;
2850 fitAC[6]=chi2AC;
2851
2852 if (fStreamLevel>2){
2853 //laser track information
2854 (*GetDebugStreamer()) << "DriftV" <<
2855 "iter=" << i <<
2856 "run=" << fRunNumber <<
2857 "time=" << timestamp <<
2858 "fitA.=" << &fitA <<
2859 "fitC.=" << &fitC <<
2860 "fitAC.=" << &fitAC <<
2861 "\n";
2862
880c3382 2863
ef7f7670 2864 }
880c3382 2865
2866 }
78f17711 2867//-----
2868
2869
2870 //Parameters: 0 linear offset (global)
2871 // 1 mean drift velocity correction factor
2872 // 2 relative global y gradient
2873 // 3 radial deformation
2874 // 4 IROC/OROC offset
2875 // 5 linear offset relative A-C
2876
2877 //get graphs
2878 TGraphErrors *grA[7];
2879 TGraphErrors *grC[7];
2880 TGraphErrors *grAC[8];
2881 TString names("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_");
2882 TString namesAC("GRAPH_MEAN_DELAY_LASER_ALL_;GRAPH_MEAN_DRIFT_LASER_ALL_;GRAPH_MEAN_GLOBALYGRADIENT_LASER_ALL_;GRAPH_MEAN_RGRADIENT_LASER_ALL_;GRAPH_MEAN_IROCOROCOFFSET_LASER_ALL_;GRAPH_MEAN_NPOINTS_LASER_ALL_;GRAPH_MEAN_CHI2_LASER_ALL_;GRAPH_MEAN_DELAYC_LASER_ALL_");
2883
2884 TObjArray *arrNames=names.Tokenize(";");
2885 TObjArray *arrNamesAC=namesAC.Tokenize(";");
2886
2887 //A-Side graphs
2888 for (Int_t i=0; i<7; ++i){
2889 TString grName=arrNames->UncheckedAt(i)->GetName();
2890 grName+="A";
2891 grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2892 if (!grA[i]){
2893 grA[i]=new TGraphErrors;
2894 grA[i]->SetName(grName.Data());
2895 grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2896 fArrFitGraphs->Add(grA[i]);
2897 }
2898// Int_t ipoint=grA[i]->GetN();
2899 Int_t ipoint=burst;
2900 grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2901 grA[i]->SetPointError(ipoint,60,.0001);
2902 if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2903 }
2904
2905 //C-Side graphs
2906 for (Int_t i=0; i<7; ++i){
2907 TString grName=arrNames->UncheckedAt(i)->GetName();
2908 grName+="C";
2909 grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2910 if (!grC[i]){
2911 grC[i]=new TGraphErrors;
2912 grC[i]->SetName(grName.Data());
2913 grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2914 fArrFitGraphs->Add(grC[i]);
2915 }
2916// Int_t ipoint=grC[i]->GetN();
2917 Int_t ipoint=burst;
2918 grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2919 grC[i]->SetPointError(ipoint,60,.0001);
2920 if (i<4) grC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2921 }
2922
2923 //AC-Side graphs
2924 for (Int_t i=0; i<8; ++i){
2925 TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2926 grName+="AC";
2927 grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2928 if (!grAC[i]){
2929 grAC[i]=new TGraphErrors;
2930 grAC[i]->SetName(grName.Data());
2931 grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2932 fArrFitGraphs->Add(grAC[i]);
2933 }
2934// Int_t ipoint=grAC[i]->GetN();
2935 Int_t ipoint=burst;
2936 grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2937 grAC[i]->SetPointError(ipoint,60,.0001);
2938 if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2939 }
2940
2941 if (fDebugLevel>10){
2942 printf("A side fit parameters:\n");
2943 fitA.Print();
2944 printf("\nC side fit parameters:\n");
2945 fitC.Print();
2946 printf("\nAC side fit parameters:\n");
2947 fitAC.Print();
2948 }
2949 delete arrNames;
2950 delete arrNamesAC;
2951}
2952
2953//_____________________________________________________________________
2954Double_t AliTPCCalibCE::SetBurstHnDrift()
2955{
2956 //
2957 // Create a new THnSparse for the current burst
2958 // return the time of the current burst
2959 //
2960 Int_t i=0;
2961 for(i=0; i<fTimeBursts.GetNrows(); ++i){
2962 if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2963 if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2964 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2965 return fTimeBursts(i);
2966 }
2967 }
2968
2969 CreateDVhist();
2970 fArrHnDrift.AddAt(fHnDrift,i);
2971 fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
2972// fPeaks[4]=0;
2973 return fTimeStamp;
2974}
2975
2976//_____________________________________________________________________
2977void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
2978{
2979 //
2980 // Write class to file
2981 // option can be specified in the dir option:
2982 // options:
2983 // name=<objname>: the name of the calibration object in file will be <objname>
2984 // type=<type>: the saving type:
2985 // 0 - write the complte object
2986 // 1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
2987 // 2 - like 2, but in addition delete objects that will most probably not be used for calibration
2988 // 3 - store only calibration output, don't store the reference histograms
2989 // and THnSparse (requires Analyse called before)
2990 //
2991 // NOTE: to read the object back, the ReadFromFile function should be used
2992 //
2993
2994 TString sDir(dir);
2995 TString objName=GetName();
2996 Int_t type=0;
2997
2998 //get options
2999 TObjArray *arr=sDir.Tokenize(",");
3000 TIter next(arr);
3001 TObjString *s;
3002 while ( (s=(TObjString*)next()) ){
3003 TString optString=s->GetString();
3004 optString.Remove(TString::kBoth,' ');
3005 if (optString.BeginsWith("name=")){
3006 objName=optString.ReplaceAll("name=","");
3007 }
3008 if (optString.BeginsWith("type=")){
3009 optString.ReplaceAll("type=","");
3010 type=optString.Atoi();
3011 }
3012 }
3013
3014 if (type==1||type==2) {
3015 //delete most probably not needed stuff
3016 //This requires Analyse to be called after reading the object from file
3017 fCalRocArrayT0.Delete();
3018 fCalRocArrayT0Err.Delete();
3019 fCalRocArrayQ.Delete();
3020 fCalRocArrayRMS.Delete();
3021 fCalRocArrayOutliers.Delete();
3022 }
3023 if (type==2||type==3){
3024 fParamArrayEventPol1.Delete();
3025 fParamArrayEventPol2.Delete();
3026 }
3027
3028 TObjArray histoQArray(72);
3029 TObjArray histoT0Array(72);
3030 TObjArray histoRMSArray(72);
3031 TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3032
3033 //save all large 2D histograms in separte pointers
3034 //to have a smaller memory print when saving the object
3035 if (type==1||type==2||type==3){
3036 for (Int_t i=0; i<72; ++i){
3037 histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3038 histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3039 histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3040 }
3041 fHistoQArray.SetOwner(kFALSE);
3042 fHistoT0Array.SetOwner(kFALSE);
3043 fHistoRMSArray.SetOwner(kFALSE);
3044 fHistoQArray.Clear();
3045 fHistoT0Array.Clear();
3046 fHistoRMSArray.Clear();
3047
3048 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3049 arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3050 }
3051 fArrHnDrift.SetOwner(kFALSE);
3052 fArrHnDrift.Clear();
3053 }
3054
3055
3056 TDirectory *backup = gDirectory;
3057
3058 TFile f(filename,"recreate");
3059 this->Write(objName.Data());
3060 if (type==1||type==2) {
3061 histoQArray.Write("histoQArray",TObject::kSingleKey);
3062 histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3063 histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3064 arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3065 }
3066
3067 f.Save();
3068 f.Close();
3069
3070 //move histograms back to the object
3071 if (type==1||type==2){
3072 for (Int_t i=0; i<72; ++i){
3073 fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3074 fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3075 fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3076 }
3077 fHistoQArray.SetOwner(kTRUE);
3078 fHistoT0Array.SetOwner(kTRUE);
3079 fHistoRMSArray.SetOwner(kTRUE);
3080
3081 for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3082 fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3083 }
3084 fArrHnDrift.SetOwner(kTRUE);
3085 }
3086
3087 if ( backup ) backup->cd();
3088}
3089//_____________________________________________________________________
3090AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3091{
3092 //
3093 // Read object from file
3094 // Handle properly if the histogram arrays were stored separately
3095 // call Analyse to make sure to have the calibration relevant information in the object
3096 //
3097
3098 TFile f(filename);
3099 if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3100 TList *l=f.GetListOfKeys();
3101 TIter next(l);
3102 TKey *key=0x0;
3103 TObject *o=0x0;
3104
3105 AliTPCCalibCE *ce=0x0;
3106 TObjArray *histoQArray=0x0;
3107 TObjArray *histoT0Array=0x0;
3108 TObjArray *histoRMSArray=0x0;
3109 TObjArray *arrHnDrift=0x0;
3110
3111 while ( (key=(TKey*)next()) ){
3112 o=key->ReadObj();
3113 if ( o->IsA()==AliTPCCalibCE::Class() ){
3114 ce=(AliTPCCalibCE*)o;
3115 } else if ( o->IsA()==TObjArray::Class() ){
3116 TString name=key->GetName();
3117 if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3118 if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3119 if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3120 if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3121 }
3122 }
3123
3124 if (ce){
3125 //move histograms back to the object
3126 TH1* hist=0x0;
3127 if (histoQArray){
3128 for (Int_t i=0; i<72; ++i){
3129 hist=(TH1*)histoQArray->UncheckedAt(i);
3130 if (hist) hist->SetDirectory(0x0);
3131 ce->fHistoQArray.AddAt(hist,i);
3132 }
3133 ce->fHistoQArray.SetOwner(kTRUE);
3134 }
3135
3136 if (histoT0Array) {
3137 for (Int_t i=0; i<72; ++i){
3138 hist=(TH1*)histoT0Array->UncheckedAt(i);
3139 if (hist) hist->SetDirectory(0x0);
3140 ce->fHistoT0Array.AddAt(hist,i);
3141 }
3142 ce->fHistoT0Array.SetOwner(kTRUE);
3143 }
3144
3145 if (histoRMSArray){
3146 for (Int_t i=0; i<72; ++i){
3147 hist=(TH1*)histoRMSArray->UncheckedAt(i);
3148 if (hist) hist->SetDirectory(0x0);
3149 ce->fHistoRMSArray.AddAt(hist,i);
3150 }
3151 ce->fHistoRMSArray.SetOwner(kTRUE);
3152 }
3153
3154 if (arrHnDrift){
3155 for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3156 THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3157 ce->fArrHnDrift.AddAt(hSparse,i);
3158 }
3159 }
3160
3161 ce->Analyse();
880c3382 3162 }
78f17711 3163 f.Close();
3164
3165 return ce;
75d8233f 3166}