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