1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 ////////////////////////////////////////////////////////////////////////////////////////
20 // Implementation of the TPC Central Electrode calibration //
22 // Origin: Jens Wiechula, Marian Ivanov J.Wiechula@gsi.de, Marian.Ivanov@cern.ch //
24 ////////////////////////////////////////////////////////////////////////////////////////
27 // *************************************************************************************
28 // * Class Description *
29 // *************************************************************************************
32 <h4>The AliTPCCalibCE class is used to get calibration data from the Central Electrode
33 using laser runs.</h4>
35 The information retrieved is
36 <ul style="list-style-type: square;">
37 <li>Time arrival from the CE</li>
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>
49 <h3><a name="working">I. Working principle</a></h3>
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>
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
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>
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>
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>
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>
90 <h4>At the end of each event the EndEvent() function is called</h4>
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>
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>
114 <h3><a name="user">II. User interface for filling data</a></h3>
116 <h4>To Fill information one of the following functions can be used:</h4>
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>
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>
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>
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>
148 <h4>It is also possible to merge two independently taken calibrations using the function</h4>
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>
160 <h4>example: filling data using root raw data:</h4>
162 void fillCE(Char_t *filename)
164 rawReader = new AliRawReaderRoot(fileName);
165 if ( !rawReader ) return;
166 AliTPCCalibCE *calib = new AliTPCCalibCE;
167 while (rawReader->NextEvent()){
168 calib->ProcessEvent(rawReader);
171 calib->DumpToFile("CEData.root");
177 <h3><a name="info">III. What kind of information is stored and how to retrieve it</a></h4>
179 <h4><a name="info:stored">III.1 Stored information</a></h4>
180 <ul style="list-style-type: none;">
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>
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>
198 <li>For each event the following information is stored:</li>
200 <ul style="list-style-type: square;">
201 <li>event time ( TVectorD fVEventTime )</li>
202 <li>event id ( TVectorD fVEventNumber )</li>
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>
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>
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>
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>
237 <h4>example for visualisation:</h4>
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");
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
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");
258 //////////////////////////////////////////////////////////////////////////////////////
262 #include <TObjArray.h>
268 #include <TVectorF.h>
269 #include <TVectorD.h>
270 #include <TVector3.h>
271 #include <TMatrixD.h>
274 #include <TGraphErrors.h>
277 #include <TDirectory.h>
280 #include <TCollection.h>
281 #include <TTimeStamp.h>
287 #include "AliRawReader.h"
288 #include "AliRawReaderRoot.h"
289 #include "AliRawReaderDate.h"
290 #include "AliRawEventHeaderBase.h"
291 #include "AliTPCRawStream.h"
292 #include "AliTPCRawStreamFast.h"
293 #include "AliTPCCalROC.h"
294 #include "AliTPCCalPad.h"
295 #include "AliTPCROC.h"
296 #include "AliTPCParam.h"
297 #include "AliTPCCalibCE.h"
298 #include "AliMathBase.h"
299 #include "AliTPCTransform.h"
300 #include "AliTPCLaserTrack.h"
301 #include "TTreeStream.h"
303 #include "AliCDBManager.h"
304 #include "AliCDBEntry.h"
307 ClassImp(AliTPCCalibCE)
310 AliTPCCalibCE::AliTPCCalibCE() :
311 AliTPCCalibRawBase(),
325 fNoiseThresholdMax(5.),
326 fNoiseThresholdSum(8.),
327 fIsZeroSuppressed(kFALSE),
330 fParam(new AliTPCParam),
336 fCalRocArrayT0Err(72),
339 fCalRocArrayOutliers(72),
347 fParamArrayEventPol1(72),
348 fParamArrayEventPol2(72),
349 fTMeanArrayEvent(72),
350 fQMeanArrayEvent(72),
357 fPadTimesArrayEvent(72),
359 fPadRMSArrayEvent(72),
360 fPadPedestalArrayEvent(72),
370 fVTime0OffsetCounter(72),
373 fCurrentCETimeRef(0),
383 // AliTPCSignal default constructor
385 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
389 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
390 for (Int_t i=0;i<5;++i){
394 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
396 //_____________________________________________________________________
397 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
398 AliTPCCalibRawBase(sig),
399 fNbinsT0(sig.fNbinsT0),
400 fXminT0(sig.fXminT0),
401 fXmaxT0(sig.fXmaxT0),
402 fNbinsQ(sig.fNbinsQ),
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),
417 fParam(new AliTPCParam),
423 fCalRocArrayT0Err(72),
426 fCalRocArrayOutliers(72),
430 fMeanT0rms(sig.fMeanT0rms),
431 fMeanQrms(sig.fMeanQrms),
432 fMeanRMSrms(sig.fMeanRMSrms),
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),
444 fPadTimesArrayEvent(72),
446 fPadRMSArrayEvent(72),
447 fPadPedestalArrayEvent(72),
457 fVTime0OffsetCounter(72),
460 fCurrentCETimeRef(0),
461 fProcessOld(sig.fProcessOld),
462 fProcessNew(sig.fProcessNew),
463 fAnalyseNew(sig.fAnalyseNew),
470 // AliTPCSignal copy constructor
472 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
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);
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);
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);
490 TH2S *hNew = new TH2S(*hQ);
491 hNew->SetDirectory(0);
492 fHistoQArray.AddAt(hNew,iSec);
495 TH2S *hNew = new TH2S(*hT0);
496 hNew->SetDirectory(0);
497 fHistoT0Array.AddAt(hNew,iSec);
500 TH2S *hNew = new TH2S(*hRMS);
501 hNew->SetDirectory(0);
502 fHistoRMSArray.AddAt(hNew,iSec);
506 //copy fit parameters event by event
508 for (Int_t iSec=0; iSec<72; ++iSec){
509 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
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);
518 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
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);
527 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
528 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
530 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
532 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
536 fVEventTime.ResizeTo(sig.fVEventTime);
537 fVEventNumber.ResizeTo(sig.fVEventNumber);
538 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
539 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
543 for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
544 TObject *o=sig.fArrHnDrift.UncheckedAt(i);
546 TObject *newo=o->Clone("fHnDrift");
547 fArrHnDrift.AddAt(newo,i);
548 if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
552 for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
553 fTimeBursts[i]=sig.fTimeBursts[i];
556 for (Int_t i=0;i<5;++i){
557 fPeaks[i]=sig.fPeaks[i];
558 fPeakWidths[i]=sig.fPeakWidths[i];
560 if (sig.fArrFitGraphs) {
561 fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
562 fArrFitGraphs->SetOwner();
565 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
568 //_____________________________________________________________________
569 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
570 AliTPCCalibRawBase(),
584 fNoiseThresholdMax(5.),
585 fNoiseThresholdSum(8.),
586 fIsZeroSuppressed(kFALSE),
589 fParam(new AliTPCParam),
595 fCalRocArrayT0Err(72),
598 fCalRocArrayOutliers(72),
606 fParamArrayEventPol1(72),
607 fParamArrayEventPol2(72),
608 fTMeanArrayEvent(72),
609 fQMeanArrayEvent(72),
616 fPadTimesArrayEvent(72),
618 fPadRMSArrayEvent(72),
619 fPadPedestalArrayEvent(72),
629 fVTime0OffsetCounter(72),
632 fCurrentCETimeRef(0),
642 // constructor which uses a tmap as input to set some specific parameters
644 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
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();
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();
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();
665 if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
666 if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
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();
672 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
673 for (Int_t i=0;i<5;++i){
679 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
682 //_____________________________________________________________________
683 AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
686 // assignment operator
688 if (&source == this) return *this;
689 new (this) AliTPCCalibCE(source);
693 //_____________________________________________________________________
694 AliTPCCalibCE::~AliTPCCalibCE()
700 fCalRocArrayT0.Delete();
701 fCalRocArrayT0Err.Delete();
702 fCalRocArrayQ.Delete();
703 fCalRocArrayRMS.Delete();
704 fCalRocArrayOutliers.Delete();
706 fHistoQArray.Delete();
707 fHistoT0Array.Delete();
708 fHistoRMSArray.Delete();
710 fHistoTmean.Delete();
712 fParamArrayEventPol1.Delete();
713 fParamArrayEventPol2.Delete();
714 fTMeanArrayEvent.Delete();
715 fQMeanArrayEvent.Delete();
717 fPadTimesArrayEvent.Delete();
718 fPadQArrayEvent.Delete();
719 fPadRMSArrayEvent.Delete();
720 fPadPedestalArrayEvent.Delete();
722 fArrHnDrift.SetOwner();
723 fArrHnDrift.Delete();
726 fArrFitGraphs->SetOwner();
727 delete fArrFitGraphs;
730 //_____________________________________________________________________
731 Int_t AliTPCCalibCE::Update(const Int_t icsector,
734 const Int_t icTimeBin,
735 const Float_t csignal)
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
743 if (!fProcessOld) return 0;
746 if (icRow<0) return 0;
747 if (icPad<0) return 0;
748 if (icTimeBin<0) return 0;
749 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
751 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
753 //init first pad and sector in this event
754 if ( fCurrentChannel == -1 ) {
756 fCurrentChannel = iChannel;
757 fCurrentSector = icsector;
761 //process last pad if we change to a new one
762 if ( iChannel != fCurrentChannel ){
764 fLastSector=fCurrentSector;
765 fCurrentChannel = iChannel;
766 fCurrentSector = icsector;
770 //fill signals for current pad
771 fPadSignal[icTimeBin]=csignal;
772 if ( csignal > fMaxPadSignal ){
773 fMaxPadSignal = csignal;
774 fMaxTimeBin = icTimeBin;
779 //_____________________________________________________________________
780 void 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)
784 // new filling method to fill the THnSparse histogram
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;
794 UShort_t timeBin = (UShort_t)startTimeBin;
795 //skip first laser layer
796 if (timeBin<200) return;
797 Double_t timeBurst=SetBurstHnDrift();
799 //after 1 event setup peak ranges
800 if (fNevents==1&&fPeaks[4]==0) {
802 fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
805 fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
808 // After the first event only fill every 5th bin in a row with the CE information
810 if (fPeaks[4]>100&&TMath::Abs((Short_t)fPeaks[4]-(Short_t)timeBin)<(Short_t)fPeakWidths[4]){
817 if (!IsPeakInRange(timeBin+length/2)) return;
819 Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
820 (Double_t)padFill,(Double_t)timeBin,timeBurst};
822 for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
823 Float_t sig=(Float_t)signal[iTimeBin];
825 fHnDrift->Fill(x,sig);
829 //_____________________________________________________________________
830 void AliTPCCalibCE::FindLaserLayers()
833 // Find the laser layer positoins
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);
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.);
849 Int_t timepos=fPeaks[4]-2*fPeakWidths[4];
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);
857 timepos=fPeaks[i]-width/2;
860 //check width and reset peak if >100
861 for (Int_t i=0; i<5; ++i){
862 if (fPeakWidths[i]>100) {
871 //_____________________________________________________________________
872 void AliTPCCalibCE::FindPedestal(Float_t part)
875 // find pedestal and noise for the current pad. Use either database or
876 // truncated mean with part*100%
878 Bool_t noPedestal = kTRUE;
880 //use pedestal database if set
881 if (fPedestalTPC&&fPadNoiseTPC){
882 //only load new pedestals if the sector has changed
883 if ( fCurrentSector!=fLastSector ){
884 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
885 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
888 if ( fPedestalROC&&fPadNoiseROC ){
889 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
890 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
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
901 if ( fIsZeroSuppressed ) return;
902 const Int_t kPedMax = 100; //maximum pedestal value
911 UShort_t histo[kPedMax];
912 memset(histo,0,kPedMax*sizeof(UShort_t));
914 //fill pedestal histogram
915 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
916 padSignal = fPadSignal[i];
917 if (padSignal<=0) continue;
918 if (padSignal>max && i>10) {
922 if (padSignal>kPedMax-1) continue;
923 histo[int(padSignal+0.5)]++;
927 for (Int_t i=1; i<kPedMax; ++i){
928 if (count1<count0*0.5) median=i;
933 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
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);
949 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
955 //_____________________________________________________________________
956 void AliTPCCalibCE::UpdateCETimeRef()
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;
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));
974 //_____________________________________________________________________
975 void AliTPCCalibCE::FindCESignal(TVectorD ¶m, Float_t &qSum, const TVectorF maxima)
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
983 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
985 Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise; // threshold for the signal sum
986 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
987 const Int_t kCemax = fPeakIntPlus;
989 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
991 // find maximum closest to the sector mean from the last event
992 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
993 // get sector mean of last event
994 Float_t tmean = fCurrentCETimeRef;
995 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
996 minDist = tmean-maxima[imax];
997 cemaxpos = (Int_t)maxima[imax];
1000 // printf("L1 phase TB: %f\n",GetL1PhaseTB());
1002 ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
1003 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
1004 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
1005 Float_t signal = fPadSignal[i]-fPadPedestal;
1007 ceTime+=signal*(i+0.5);
1008 ceRMS +=signal*(i+0.5)*(i+0.5);
1014 if (ceQmax&&ceQsum>ceSumThreshold) {
1016 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
1017 ceTime-=GetL1PhaseTB();
1018 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
1019 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
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.;
1030 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
1031 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
1043 //_____________________________________________________________________
1044 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
1047 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
1048 // and 'tplus' timebins after 'pos'
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
1056 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
1060 //_____________________________________________________________________
1061 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
1064 // Find local maxima on the pad signal and Histogram them
1066 Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal
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) ){
1072 if (count<maxima.GetNrows()){
1073 maxima.GetMatrixArray()[count++]=i;
1074 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
1075 i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin fPeakDetMinus+fPeakDetPlus-1
1080 //_____________________________________________________________________
1081 void AliTPCCalibCE::ProcessPad()
1084 // Process data of current pad
1088 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
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
1091 FindLocalMaxima(maxima);
1092 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
1094 UpdateCETimeRef(); // update the time refenrence for the current sector
1095 if ( fCurrentCETimeRef<1e-30 ) return; //return if we don't have time 0 info, eg if only one side has laser
1098 FindCESignal(param, qSum, maxima);
1100 Double_t meanT = param[1];
1101 Double_t sigmaT = param[2];
1103 //Fill Event T0 counter
1104 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
1107 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
1109 //Fill RMS histogram
1110 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
1113 //Fill debugging info
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;
1122 //_____________________________________________________________________
1123 void AliTPCCalibCE::EndEvent()
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
1131 if (fProcessNew) ++fNevents;
1135 //check if last pad has allready been processed, if not do so
1136 if ( fMaxTimeBin>-1 ) ProcessPad();
1138 AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
1141 TMatrixD dummy(3,3);
1142 // TVectorF vMeanTime(72);
1143 // TVectorF vMeanQ(72);
1144 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
1145 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
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];
1156 if ( time0SideCount[0] >0 )
1157 time0Side[0]/=time0SideCount[0];
1158 if ( time0SideCount[1] >0 )
1159 time0Side[1]/=time0SideCount[1];
1160 // end find time0 offset
1161 AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
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 ){
1172 AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1176 Double_t entries = hMeanT->GetEffectiveEntries();
1178 Short_t *arr = hMeanT->GetArray()+1;
1180 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1182 if ( sum>=(entries/2.) ) break;
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);
1191 // check boundaries for ebye info of mean time
1192 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1193 Int_t vSize=vMeanTime->GetNrows();
1194 if ( vSize < fNevents+1 ){
1195 vMeanTime->ResizeTo(vSize+100);
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);
1204 fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1205 fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
1207 vMeanTime->GetMatrixArray()[fNevents]=median;
1211 TVectorF *vTimes = GetPadTimesEvent(iSec);
1212 if ( !vTimes ) continue; //continue if no time information for this sector is available
1214 AliTPCCalROC calIrocOutliers(0);
1215 AliTPCCalROC calOrocOutliers(36);
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);
1224 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1225 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1227 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1228 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
1230 //set values for temporary roc calibration class
1232 calIroc->SetValue(iChannel, time);
1233 if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
1236 calOroc->SetValue(iChannel, time);
1237 if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
1240 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
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 );
1247 //------------------------------- Debug start ------------------------------
1248 if ( GetStreamLevel()>0 ){
1249 TTreeSRedirector *streamer=GetDebugStreamer();
1255 Float_t q = (*GetPadQEvent(iSec))[iChannel];
1256 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1258 UInt_t channel=iChannel;
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);
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);
1271 // for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1272 // h1->Fill(i,fPadSignal(i));
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" <<
1279 "Event=" << fNevents <<
1280 "RunNumber=" << fRunNumber <<
1281 "TimeStamp=" << fTimeStamp <<
1282 "Sector="<< sector <<
1286 "PadSec="<< channel <<
1287 "Time0Sec=" << t0Sec <<
1288 "Time0Side=" << t0Side <<
1292 "MeanQ=" << meanQ <<
1293 // "hist.=" << h1 <<
1299 //----------------------------- Debug end ------------------------------
1300 }// end channel loop
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);
1312 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1314 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1315 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1317 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1318 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1321 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1322 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1325 //------------------------------- Debug start ------------------------------
1326 if ( GetStreamLevel()>0 ){
1327 TTreeSRedirector *streamer=GetDebugStreamer();
1329 (*streamer) << "DataRoc" <<
1330 // "Event=" << fEvent <<
1331 "RunNumber=" << fRunNumber <<
1332 "TimeStamp=" << fTimeStamp <<
1334 "hMeanT.=" << hMeanT <<
1335 "median=" << median <<
1336 "paramPol1.=" << ¶mPol1 <<
1337 "paramPol2.=" << ¶mPol2 <<
1338 "matPol1.=" << &matPol1 <<
1339 "matPol2.=" << &matPol2 <<
1340 "chi2Pol1=" << chi2Pol1 <<
1341 "chi2Pol2=" << chi2Pol2 <<
1346 //------------------------------- Debug end ------------------------------
1349 //return if no sector has a valid mean time
1350 if ( nSecMeanT == 0 ) return;
1353 // fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1354 // fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1355 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1356 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1357 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1359 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1360 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1363 fOldRunNumber = fRunNumber;
1367 AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1369 //_____________________________________________________________________
1370 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1371 Int_t nbinsY, Float_t ymin, Float_t ymax,
1372 const Char_t *type, Bool_t force)
1375 // return pointer to TH2S histogram of 'type'
1376 // if force is true create a new histogram if it doesn't exist allready
1378 if ( !force || arr->UncheckedAt(sector) )
1379 return (TH2S*)arr->UncheckedAt(sector);
1381 // if we are forced and histogram doesn't exist yet create it
1382 // new histogram with Q calib information. One value for each pad!
1383 TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1385 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1386 hist->SetDirectory(0);
1387 arr->AddAt(hist,sector);
1390 //_____________________________________________________________________
1391 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1394 // return pointer to T0 histogram
1395 // if force is true create a new histogram if it doesn't exist allready
1397 TObjArray *arr = &fHistoT0Array;
1398 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1400 //_____________________________________________________________________
1401 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1404 // return pointer to Q histogram
1405 // if force is true create a new histogram if it doesn't exist allready
1407 TObjArray *arr = &fHistoQArray;
1408 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1410 //_____________________________________________________________________
1411 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1414 // return pointer to Q histogram
1415 // if force is true create a new histogram if it doesn't exist allready
1417 TObjArray *arr = &fHistoRMSArray;
1418 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1420 //_____________________________________________________________________
1421 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1422 const Char_t *type, Bool_t force)
1425 // return pointer to TH1S histogram
1426 // if force is true create a new histogram if it doesn't exist allready
1428 if ( !force || arr->UncheckedAt(sector) )
1429 return (TH1S*)arr->UncheckedAt(sector);
1431 // if we are forced and histogram doesn't yes exist create it
1432 // new histogram with calib information. One value for each pad!
1433 TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1434 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1435 hist->SetDirectory(0);
1436 arr->AddAt(hist,sector);
1439 //_____________________________________________________________________
1440 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1443 // return pointer to Q histogram
1444 // if force is true create a new histogram if it doesn't exist allready
1446 TObjArray *arr = &fHistoTmean;
1447 return GetHisto(sector, arr, "LastTmean", force);
1449 //_____________________________________________________________________
1450 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1453 // return pointer to Pad Info from 'arr' for the current event and sector
1454 // if force is true create it if it doesn't exist allready
1456 if ( !force || arr->UncheckedAt(sector) )
1457 return (TVectorF*)arr->UncheckedAt(sector);
1459 TVectorF *vect = new TVectorF(size);
1460 arr->AddAt(vect,sector);
1463 //_____________________________________________________________________
1464 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1467 // return pointer to Pad Times Array for the current event and sector
1468 // if force is true create it if it doesn't exist allready
1470 TObjArray *arr = &fPadTimesArrayEvent;
1471 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1473 //_____________________________________________________________________
1474 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1477 // return pointer to Pad Q Array for the current event and sector
1478 // if force is true create it if it doesn't exist allready
1479 // for debugging purposes only
1482 TObjArray *arr = &fPadQArrayEvent;
1483 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1485 //_____________________________________________________________________
1486 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1489 // return pointer to Pad RMS Array for the current event and sector
1490 // if force is true create it if it doesn't exist allready
1491 // for debugging purposes only
1493 TObjArray *arr = &fPadRMSArrayEvent;
1494 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1496 //_____________________________________________________________________
1497 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1500 // return pointer to Pad RMS Array for the current event and sector
1501 // if force is true create it if it doesn't exist allready
1502 // for debugging purposes only
1504 TObjArray *arr = &fPadPedestalArrayEvent;
1505 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1507 //_____________________________________________________________________
1508 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1511 // return pointer to the EbyE info of the mean arrival time for 'sector'
1512 // if force is true create it if it doesn't exist allready
1514 TObjArray *arr = &fTMeanArrayEvent;
1515 return GetVectSector(sector,arr,100,force);
1517 //_____________________________________________________________________
1518 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
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
1524 TObjArray *arr = &fQMeanArrayEvent;
1525 return GetVectSector(sector,arr,100,force);
1527 //_____________________________________________________________________
1528 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1531 // return pointer to ROC Calibration
1532 // if force is true create a new histogram if it doesn't exist allready
1534 if ( !force || arr->UncheckedAt(sector) )
1535 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1537 // if we are forced and histogram doesn't yes exist create it
1539 // new AliTPCCalROC for T0 information. One value for each pad!
1540 AliTPCCalROC *croc = new AliTPCCalROC(sector);
1541 arr->AddAt(croc,sector);
1544 //_____________________________________________________________________
1545 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1548 // return pointer to Time 0 ROC Calibration
1549 // if force is true create a new histogram if it doesn't exist allready
1551 TObjArray *arr = &fCalRocArrayT0;
1552 return GetCalRoc(sector, arr, force);
1554 //_____________________________________________________________________
1555 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1558 // return pointer to the error of Time 0 ROC Calibration
1559 // if force is true create a new histogram if it doesn't exist allready
1561 TObjArray *arr = &fCalRocArrayT0Err;
1562 return GetCalRoc(sector, arr, force);
1564 //_____________________________________________________________________
1565 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1568 // return pointer to T0 ROC Calibration
1569 // if force is true create a new histogram if it doesn't exist allready
1571 TObjArray *arr = &fCalRocArrayQ;
1572 return GetCalRoc(sector, arr, force);
1574 //_____________________________________________________________________
1575 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1578 // return pointer to signal width ROC Calibration
1579 // if force is true create a new histogram if it doesn't exist allready
1581 TObjArray *arr = &fCalRocArrayRMS;
1582 return GetCalRoc(sector, arr, force);
1584 //_____________________________________________________________________
1585 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1588 // return pointer to Outliers
1589 // if force is true create a new histogram if it doesn't exist allready
1591 TObjArray *arr = &fCalRocArrayOutliers;
1592 return GetCalRoc(sector, arr, force);
1594 //_____________________________________________________________________
1595 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1598 // return pointer to TObjArray of fit parameters
1599 // if force is true create a new histogram if it doesn't exist allready
1601 if ( !force || arr->UncheckedAt(sector) )
1602 return (TObjArray*)arr->UncheckedAt(sector);
1604 // if we are forced and array doesn't yes exist create it
1606 // new TObjArray for parameters
1607 TObjArray *newArr = new TObjArray;
1608 arr->AddAt(newArr,sector);
1611 //_____________________________________________________________________
1612 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1615 // return pointer to TObjArray of fit parameters from plane fit
1616 // if force is true create a new histogram if it doesn't exist allready
1618 TObjArray *arr = &fParamArrayEventPol1;
1619 return GetParamArray(sector, arr, force);
1621 //_____________________________________________________________________
1622 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1625 // return pointer to TObjArray of fit parameters from parabola fit
1626 // if force is true create a new histogram if it doesn't exist allready
1628 TObjArray *arr = &fParamArrayEventPol2;
1629 return GetParamArray(sector, arr, force);
1632 //_____________________________________________________________________
1633 void AliTPCCalibCE::CreateDVhist()
1636 // Setup the THnSparse for the drift velocity determination
1640 //roc, row, pad, timebin, timestamp
1641 TTimeStamp begin(2010,01,01,0,0,0);
1642 TTimeStamp end(2035,01,01,0,0,0);
1643 Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1645 Int_t bins[kHnBinsDV] = { 72, 96, 140, 1030, nbinsTime};
1646 Double_t xmin[kHnBinsDV] = { 0., 0., 0., 0., (Double_t)begin.GetSec()};
1647 Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1649 fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1650 fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1651 fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1652 fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1653 fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1654 fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
1658 //_____________________________________________________________________
1659 void AliTPCCalibCE::ResetEvent()
1662 // Reset global counters -- Should be called before each event is processed
1671 fPadTimesArrayEvent.Delete();
1672 fPadQArrayEvent.Delete();
1673 fPadRMSArrayEvent.Delete();
1674 fPadPedestalArrayEvent.Delete();
1676 for ( Int_t i=0; i<72; ++i ){
1677 fVTime0Offset.GetMatrixArray()[i]=0;
1678 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1679 fVMeanQ.GetMatrixArray()[i]=0;
1680 fVMeanQCounter.GetMatrixArray()[i]=0;
1683 //_____________________________________________________________________
1684 void AliTPCCalibCE::ResetPad()
1687 // Reset pad infos -- Should be called after a pad has been processed
1689 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1696 //_____________________________________________________________________
1697 void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
1700 // Merge ce to the current AliTPCCalibCE
1703 Int_t nCEevents = ce->GetNeventsProcessed();
1705 if (fProcessOld&&ce->fProcessOld){
1707 for (Int_t iSec=0; iSec<72; ++iSec){
1708 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1709 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1710 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1714 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1715 TH2S *hRefQ = GetHistoQ(iSec);
1716 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1718 TH2S *hist = new TH2S(*hRefQmerge);
1719 hist->SetDirectory(0);
1720 fHistoQArray.AddAt(hist, iSec);
1722 hRefQmerge->SetDirectory(dir);
1725 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1726 TH2S *hRefT0 = GetHistoT0(iSec);
1727 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1729 TH2S *hist = new TH2S(*hRefT0merge);
1730 hist->SetDirectory(0);
1731 fHistoT0Array.AddAt(hist, iSec);
1733 hRefT0merge->SetDirectory(dir);
1735 if ( hRefRMSmerge ){
1736 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1737 TH2S *hRefRMS = GetHistoRMS(iSec);
1738 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1740 TH2S *hist = new TH2S(*hRefRMSmerge);
1741 hist->SetDirectory(0);
1742 fHistoRMSArray.AddAt(hist, iSec);
1744 hRefRMSmerge->SetDirectory(dir);
1749 // merge time information
1752 for (Int_t iSec=0; iSec<72; ++iSec){
1753 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1754 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1755 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1756 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1758 TObjArray *arrPol1 = 0x0;
1759 TObjArray *arrPol2 = 0x0;
1760 TVectorF *vMeanTime = 0x0;
1761 TVectorF *vMeanQ = 0x0;
1764 if ( arrPol1CE && arrPol2CE ){
1765 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1766 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1767 arrPol1->Expand(fNevents+nCEevents);
1768 arrPol2->Expand(fNevents+nCEevents);
1770 if ( vMeanTimeCE && vMeanQCE ){
1771 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1772 vMeanQ = GetQMeanEvents(iSec,kTRUE);
1773 vMeanTime->ResizeTo(fNevents+nCEevents);
1774 vMeanQ->ResizeTo(fNevents+nCEevents);
1777 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1778 if ( arrPol1CE && arrPol2CE ){
1779 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1780 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1781 if ( paramPol1 && paramPol2 ){
1782 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1783 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1786 if ( vMeanTimeCE && vMeanQCE ){
1787 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1788 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1795 const TVectorD& eventTimes = ce->fVEventTime;
1796 const TVectorD& eventIds = ce->fVEventNumber;
1797 const TVectorF& time0SideA = ce->fVTime0SideA;
1798 const TVectorF& time0SideC = ce->fVTime0SideC;
1799 fVEventTime.ResizeTo(fNevents+nCEevents);
1800 fVEventNumber.ResizeTo(fNevents+nCEevents);
1801 fVTime0SideA.ResizeTo(fNevents+nCEevents);
1802 fVTime0SideC.ResizeTo(fNevents+nCEevents);
1804 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1805 Double_t evTime = eventTimes.GetMatrixArray()[iEvent];
1806 Double_t evId = eventIds.GetMatrixArray()[iEvent];
1807 Float_t t0SideA = time0SideA.GetMatrixArray()[iEvent];
1808 Float_t t0SideC = time0SideC.GetMatrixArray()[iEvent];
1810 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1811 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1812 fVTime0SideA.GetMatrixArray()[fNevents+iEvent] = t0SideA;
1813 fVTime0SideC.GetMatrixArray()[fNevents+iEvent] = t0SideC;
1817 if (fProcessNew&&ce->fProcessNew) {
1818 if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1819 AliError("Number of bursts in the instances to merge are different. No merging done!");
1821 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1822 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1823 THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1824 if (h && hce) h->Add(hce);
1825 else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
1827 //TODO: What to do with fTimeBursts???
1831 fNevents+=nCEevents; //increase event counter
1834 //_____________________________________________________________________
1835 Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1838 // Merge all objects of this type in list
1844 AliTPCCalibCE *ce=0;
1847 while ( (o=next()) ){
1848 ce=dynamic_cast<AliTPCCalibCE*>(o);
1858 //_____________________________________________________________________
1859 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1862 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1863 // or side (-1: A-Side, -2: C-Side)
1864 // xVariable: 0-event time, 1-event id, 2-internal event counter
1865 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1866 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1867 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1868 // not used for mean time and mean Q )
1869 // for an example see class description at the beginning
1872 TVectorD *xVar = 0x0;
1873 TObjArray *aType = 0x0;
1877 if ( (sector<-2) || (sector>71) ) return 0x0; //sector outside valid range
1878 if ( (xVariable<0) || (xVariable>2) ) return 0x0; //invalid x-variable
1879 if ( (fitType<0) || (fitType>3) ) return 0x0; //invalid fit type
1880 if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1881 if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1882 if ( sector<0 && fitType!=2) return 0x0; //for side wise information only mean time is available
1886 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1887 aType = &fParamArrayEventPol1;
1888 if ( aType->At(sector)==0x0 ) return 0x0;
1890 else if ( fitType==1 ){
1891 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1892 aType = &fParamArrayEventPol2;
1893 if ( aType->At(sector)==0x0 ) return 0x0;
1897 if ( xVariable == 0 ) xVar = &fVEventTime;
1898 if ( xVariable == 1 ) xVar = &fVEventNumber;
1899 if ( xVariable == 2 ) {
1900 xVar = new TVectorD(fNevents);
1901 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1904 Double_t *x = new Double_t[fNevents];
1905 Double_t *y = new Double_t[fNevents];
1907 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1909 TObjArray *events = (TObjArray*)(aType->At(sector));
1910 if ( events->GetSize()<=ievent ) break;
1911 TVectorD *v = (TVectorD*)(events->At(ievent));
1912 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1913 } else if (fitType == 2) {
1914 Double_t xValue=(*xVar)[ievent];
1916 if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1917 else if (sector==-1) yValue=fVTime0SideA(ievent);
1918 else if (sector==-2) yValue=fVTime0SideC(ievent);
1919 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1920 }else if (fitType == 3) {
1921 Double_t xValue=(*xVar)[ievent];
1922 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1923 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1927 TGraph *gr = new TGraph(npoints);
1928 //sort xVariable increasing
1929 Int_t *sortIndex = new Int_t[npoints];
1930 TMath::Sort(npoints,x,sortIndex, kFALSE);
1931 for (Int_t i=0;i<npoints;++i){
1932 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1936 if ( xVariable == 2 ) delete xVar;
1939 delete [] sortIndex;
1942 //_____________________________________________________________________
1943 void AliTPCCalibCE::Analyse()
1946 // Calculate calibration constants
1951 TVectorD paramT0(3);
1952 TVectorD paramRMS(3);
1953 TMatrixD dummy(3,3);
1955 Float_t channelCounter=0;
1960 for (Int_t iSec=0; iSec<72; ++iSec){
1961 TH2S *hT0 = GetHistoT0(iSec);
1962 if (!hT0 ) continue;
1964 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1965 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1966 AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1967 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1968 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
1970 TH2S *hQ = GetHistoQ(iSec);
1971 TH2S *hRMS = GetHistoRMS(iSec);
1973 Short_t *arrayhQ = hQ->GetArray();
1974 Short_t *arrayhT0 = hT0->GetArray();
1975 Short_t *arrayhRMS = hRMS->GetArray();
1977 UInt_t nChannels = fROC->GetNChannels(iSec);
1985 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
1988 Float_t cogTime0 = -1000;
1989 Float_t cogQ = -1000;
1990 Float_t cogRMS = -1000;
1996 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
1997 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
1998 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
2000 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2002 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2004 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2009 //outlier specifications
2010 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2017 rocQ->SetValue(iChannel, cogQ*cogQ);
2018 rocT0->SetValue(iChannel, cogTime0);
2019 rocT0Err->SetValue(iChannel, rmsT0);
2020 rocRMS->SetValue(iChannel, cogRMS);
2021 rocOut->SetValue(iChannel, cogOut);
2025 if ( GetStreamLevel() > 0 ){
2026 TTreeSRedirector *streamer=GetDebugStreamer();
2029 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2030 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2031 padc = pad-(fROC->GetNPads(iSec,row)/2);
2033 (*streamer) << "DataEnd" <<
2034 "Sector=" << iSec <<
2038 "PadSec=" << iChannel <<
2040 "T0=" << cogTime0 <<
2050 if ( channelCounter>0 ){
2051 fMeanT0rms/=channelCounter;
2052 fMeanQrms/=channelCounter;
2053 fMeanRMSrms/=channelCounter;
2055 // if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2056 // delete fDebugStreamer;
2057 // fDebugStreamer = 0x0;
2058 fVEventTime.ResizeTo(fNevents);
2059 fVEventNumber.ResizeTo(fNevents);
2060 fVTime0SideA.ResizeTo(fNevents);
2061 fVTime0SideC.ResizeTo(fNevents);
2064 if (fProcessNew&&fAnalyseNew){
2066 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2067 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2068 h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2077 // New functions that also use the laser tracks
2082 //_____________________________________________________________________
2083 void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2086 //Find the local maximums for the projections to each axis
2089 //find laser layer positoins
2090 fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2092 THnSparse *hProj=fHnDrift;
2093 Double_t posCE[4]={0.,0.,0.,0.};
2094 Double_t widthCE[4]={0.,0.,0.,0.};
2096 // if(fPeaks[4]!=0){
2097 // find central electrode position once more, separately for IROC, OROC, A-, C-Side
2099 for (Int_t i=0; i<4; ++i){
2100 hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2101 TH1 *h=fHnDrift->Projection(3);
2102 h->GetXaxis()->SetRangeUser(fPeaks[4]-fPeakWidths[4],fPeaks[4]+fPeakWidths[4]);
2103 Int_t nbinMax=h->GetMaximumBin();
2104 Double_t maximum=h->GetMaximum();
2105 // Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2106 // if (nbinMax<700||maximum<maxExpected) continue;
2107 Double_t xbinMax=h->GetBinCenter(nbinMax);
2108 TF1 fgaus("gaus","gaus",xbinMax-10,xbinMax+10);
2109 fgaus.SetParameters(maximum,xbinMax,2);
2110 fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2111 fgaus.SetParLimits(2,0.2,4.);
2112 h->Fit(&fgaus,"RQN");
2113 // Double_t deltaX=4*fgaus.GetParameter(2);
2114 // xbinMax=fgaus.GetParameter(1);
2116 posCE[i]=fgaus.GetParameter(1);
2117 widthCE[i]=4*fgaus.GetParameter(2);
2118 hProj->GetAxis(0)->SetRangeUser(0,72);
2121 //Current drift velocity
2122 Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2123 // cout<<"vdrift="<<vdrift<<endl;
2125 AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2126 //loop over all entries in the histogram
2128 for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2129 //get entry position and content
2130 Double_t adc=hProj->GetBinContent(ichunk,coord);
2133 Int_t sector = coord[0]-1;
2134 Int_t row = coord[1]-1;
2135 Int_t pad = coord[2]-1;
2136 Int_t timeBin= coord[3]-1;
2137 Double_t time = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2138 Int_t side = (sector/18)%2;
2140 // fPeaks[4]=(UInt_t)posCE[sector/18];
2141 // fPeakWidths[4]=(UInt_t)widthCE[sector/18];
2144 if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2145 if (adc < 5 ) continue;
2146 if (IsEdgePad(sector,row,pad)) continue;
2147 // if (!IsPeakInRange(timeBin)) continue;
2148 // if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2149 // if (isCE&&(sector!=0)) continue;
2151 Int_t padmin=-2, padmax=2;
2152 Int_t timemin=-2, timemax=2;
2153 Int_t minsumperpad=2;
2154 //CE or laser tracks
2156 if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2165 // Find local maximum and cogs
2167 Bool_t isMaximum=kTRUE;
2168 Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2169 Double_t cogY=0, rmsY=0;
2172 // for position calculation use
2173 for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2175 fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2177 for(Int_t itime=timemin;itime<=timemax;++itime){
2179 Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2180 Double_t val=hProj->GetBinContent(a);
2183 tbcm +=(timeBin+itime)*val;
2184 padcm+=(pad+ipad)*val;
2187 rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2188 rmspad+=(pad+ipad)*(pad+ipad)*val;
2189 rmsY +=lxyz[1]*lxyz[1]*val;
2197 if (!isMaximum) break;
2200 if (!isMaximum||!npart) continue;
2201 if (totalmass<npart*minsumperpad) continue;
2202 if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2203 //for CE we want only one pad by construction
2213 rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2214 rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2215 rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2217 Int_t cog=TMath::Nint(padcm);
2219 // timebin --> z position
2220 Float_t zlength=fParam->GetZLength(side);
2221 // Float_t timePos=tbcm+(1000-fPeaks[4])
2222 // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2223 Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2225 // local to global transformation--> x and y positions
2227 fROC->GetPositionLocal(sector,row,pad,padlxyz);
2229 Double_t gxyz[3]={padlxyz[0],cogY,gz};
2230 Double_t lxyz[3]={padlxyz[0],cogY,gz};
2231 Double_t igxyz[3]={0,0,0};
2233 t1.RotatedGlobal2Global(sector,gxyz);
2239 //find track id and index of the position in the track (row)
2242 index=row+(sector>35)*fROC->GetNRows(0);
2243 trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2245 trackID=336+((sector/18)%2);
2246 index= fROC->GetRowIndexes(sector)[row]+pad; // global pad position in sector
2248 index+=(sector%18)*fROC->GetNChannels(sector);
2250 index+=18*fROC->GetNChannels(0);
2251 index+=(sector%18)*fROC->GetNChannels(sector);
2253 //TODO: find out about the multiple peaks in the CE
2254 // mindist=TMath::Abs(fPeaks[4]-tbcm);
2258 // fill track vectors
2260 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2261 Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2264 // travel time effect of light includes
2266 Double_t raylength=ltr->GetRayLength();
2267 Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2268 ltr->GetXYZ(globmir);
2271 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2272 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2273 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2276 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2277 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2278 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2282 if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2283 ltr->fVecSec->GetMatrixArray()[index]=sector;
2284 ltr->fVecP2->GetMatrixArray()[index]=0;
2285 ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2286 ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2287 ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2288 ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2289 ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2290 ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2291 ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2292 // ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2294 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2295 ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2296 igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2297 igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2298 igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2302 if (fStreamLevel>4){
2303 (*GetDebugStreamer()) << "clusters" <<
2304 "run=" << fRunNumber <<
2305 "timestamp=" << timestamp <<
2306 "burst=" << burst <<
2312 "timebin=" << timeBin <<
2313 "cogCE=" << posCE[sector/18] <<
2314 "withCE=" << widthCE[sector/18] <<
2315 "index=" << index <<
2317 "padcm=" << padcm <<
2318 "rmspad=" << rmspad <<
2321 "rmstb=" << rmstb <<
2325 "lx=" << padlxyz[0]<<
2327 "lypad=" << padlxyz[1]<<
2334 "igx=" << igxyz[0] <<
2335 "igy=" << igxyz[1] <<
2336 "igz=" << igxyz[2] <<
2338 "mind=" << mindist <<
2340 "trackid=" << trackID <<
2341 "trackid2=" << trackID2 <<
2342 "npart=" << npart <<
2344 } // end stream levelmgz.fElements
2350 //_____________________________________________________________________
2351 void AliTPCCalibCE::AnalyseTrack()
2354 // Analyse the tracks
2358 AliTPCLaserTrack::LoadTracks();
2359 // AliTPCParam *param=0x0;
2361 // AliCDBManager *man=AliCDBManager::Instance();
2362 // if (man->GetDefaultStorage()){
2363 // AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2365 // entry->SetOwner(kTRUE);
2366 // param = (AliTPCParam*)(entry->GetObject()->Clone());
2370 // if (fParam) delete fParam;
2373 // AliError("Could not get updated AliTPCParam from OCDB!!!");
2376 //Measured and ideal laser tracks
2377 TObjArray* arrMeasured = SetupMeasured();
2378 TObjArray* arrIdeal = AliTPCLaserTrack::GetTracks();
2379 AddCEtoIdeal(arrIdeal);
2381 //find bursts and loop over them
2382 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2383 Double_t timestamp=fTimeBursts[iburst];
2384 AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2385 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2386 if (!fHnDrift) continue;
2387 UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2388 if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2389 fBinsLastAna[iburst]=entries;
2391 for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2392 // if (iburst==0) FindLaserLayers();
2394 //reset laser tracks
2395 ResetMeasured(arrMeasured);
2397 //find clusters and associate to the tracks
2398 FindLocalMaxima(arrMeasured, timestamp, iburst);
2400 //calculate drift velocity
2401 CalculateDV(arrIdeal,arrMeasured,iburst);
2403 //Dump information to file if requested
2404 if (fStreamLevel>2){
2405 printf("make tree\n");
2406 //laser track information
2408 for (Int_t itrack=0; itrack<338; ++itrack){
2409 TObject *iltr=arrIdeal->UncheckedAt(itrack);
2410 TObject *mltr=arrMeasured->UncheckedAt(itrack);
2411 (*GetDebugStreamer()) << "tracks" <<
2412 "run=" << fRunNumber <<
2413 "time=" << timestamp <<
2414 "burst="<< iburst <<
2421 if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2424 //_____________________________________________________________________
2425 Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2426 const Double_t *peakposloc, Int_t &itrackMin2)
2429 // Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2433 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2434 TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2439 Int_t lastbeam=336/2;
2440 if ( (sector/18)%2 ) {
2447 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2448 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2449 // if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2450 vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2451 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2452 if (TMath::Abs(deltaZ)>40) continue;
2453 vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2454 vSt.RotateZ(ltr->GetAlpha());
2455 vDir.RotateZ(ltr->GetAlpha());
2457 Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2466 Float_t mindist2=10;
2467 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2468 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2469 if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2471 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2472 if (TMath::Abs(deltaZ)>40) continue;
2474 Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2475 if (dist>1) continue;
2487 //_____________________________________________________________________
2488 Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2491 // return true if pad is on the edge of a row
2494 if ( pad == edge1 ) return kTRUE;
2495 Int_t edge2 = fROC->GetNPads(sector,row)-1;
2496 if ( pad == edge2 ) return kTRUE;
2501 //_____________________________________________________________________
2502 TObjArray* AliTPCCalibCE::SetupMeasured()
2505 // setup array of measured laser tracks and CE
2508 TObjArray *arrIdeal = AliTPCLaserTrack::GetTracks();
2509 TObjArray *arrMeasured = new TObjArray(338);
2510 arrMeasured->SetOwner();
2511 for(Int_t itrack=0;itrack<336;itrack++){
2512 arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2516 AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2519 ltrce->fVecSec=new TVectorD(557568/2);
2520 ltrce->fVecP2=new TVectorD(557568/2);
2521 ltrce->fVecPhi=new TVectorD(557568/2);
2522 ltrce->fVecGX=new TVectorD(557568/2);
2523 ltrce->fVecGY=new TVectorD(557568/2);
2524 ltrce->fVecGZ=new TVectorD(557568/2);
2525 ltrce->fVecLX=new TVectorD(557568/2);
2526 ltrce->fVecLY=new TVectorD(557568/2);
2527 ltrce->fVecLZ=new TVectorD(557568/2);
2529 arrMeasured->AddAt(ltrce,336); //CE A-Side
2531 ltrce=new AliTPCLaserTrack;
2534 ltrce->fVecSec=new TVectorD(557568/2);
2535 ltrce->fVecP2=new TVectorD(557568/2);
2536 ltrce->fVecPhi=new TVectorD(557568/2);
2537 ltrce->fVecGX=new TVectorD(557568/2);
2538 ltrce->fVecGY=new TVectorD(557568/2);
2539 ltrce->fVecGZ=new TVectorD(557568/2);
2540 ltrce->fVecLX=new TVectorD(557568/2);
2541 ltrce->fVecLY=new TVectorD(557568/2);
2542 ltrce->fVecLZ=new TVectorD(557568/2);
2543 arrMeasured->AddAt(ltrce,337); //CE C-Side
2548 //_____________________________________________________________________
2549 void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2552 // reset array of measured laser tracks and CE
2554 for(Int_t itrack=0;itrack<338;itrack++){
2555 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2556 ltr->fVecSec->Zero();
2557 ltr->fVecP2->Zero();
2558 ltr->fVecPhi->Zero();
2559 ltr->fVecGX->Zero();
2560 ltr->fVecGY->Zero();
2561 ltr->fVecGZ->Zero();
2562 ltr->fVecLX->Zero();
2563 ltr->fVecLY->Zero();
2564 ltr->fVecLZ->Zero();
2568 //_____________________________________________________________________
2569 void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2572 // Add ideal CE positions to the ideal track data
2577 AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2580 ltrceA->fVecSec=new TVectorD(557568/2);
2581 ltrceA->fVecP2=new TVectorD(557568/2);
2582 ltrceA->fVecPhi=new TVectorD(557568/2);
2583 ltrceA->fVecGX=new TVectorD(557568/2);
2584 ltrceA->fVecGY=new TVectorD(557568/2);
2585 ltrceA->fVecGZ=new TVectorD(557568/2);
2586 ltrceA->fVecLX=new TVectorD(557568/2);
2587 ltrceA->fVecLY=new TVectorD(557568/2);
2588 ltrceA->fVecLZ=new TVectorD(557568/2);
2589 arr->AddAt(ltrceA,336); //CE A-Side
2591 AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2594 ltrceC->fVecSec=new TVectorD(557568/2);
2595 ltrceC->fVecP2=new TVectorD(557568/2);
2596 ltrceC->fVecPhi=new TVectorD(557568/2);
2597 ltrceC->fVecGX=new TVectorD(557568/2);
2598 ltrceC->fVecGY=new TVectorD(557568/2);
2599 ltrceC->fVecGZ=new TVectorD(557568/2);
2600 ltrceC->fVecLX=new TVectorD(557568/2);
2601 ltrceC->fVecLY=new TVectorD(557568/2);
2602 ltrceC->fVecLZ=new TVectorD(557568/2);
2603 arr->AddAt(ltrceC,337); //CE C-Side
2605 //Calculate ideal positoins
2608 Int_t channelSideA=0;
2609 Int_t channelSideC=0;
2610 Int_t channelSide=0;
2611 AliTPCLaserTrack *ltrce=0x0;
2613 for (Int_t isector=0; isector<72; ++isector){
2614 Int_t side=((isector/18)%2);
2615 for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2616 for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2617 fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2618 fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2621 channelSide=channelSideA;
2624 channelSide=channelSideC;
2627 ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2628 ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2629 ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2630 ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2631 ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2632 // ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2633 ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2634 ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2635 // ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2638 ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2639 ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2643 ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2644 ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2654 //_____________________________________________________________________
2655 void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2658 // calculate the drift velocity from the reconstructed clusters associated
2659 // to the ideal laser tracks
2660 // use two different fit scenarios: Separate fit for A- and C-Side
2661 // Common fit for A- and C-Side
2664 if (!fArrFitGraphs){
2665 fArrFitGraphs=new TObjArray;
2668 // static TLinearFitter fdriftA(5,"hyp4");
2669 // static TLinearFitter fdriftC(5,"hyp4");
2670 // static TLinearFitter fdriftAC(6,"hyp5");
2671 Double_t timestamp=fTimeBursts[burst];
2673 static TLinearFitter fdriftA(4,"hyp3");
2674 static TLinearFitter fdriftC(4,"hyp3");
2675 static TLinearFitter fdriftAC(5,"hyp4");
2676 TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2680 Float_t chi2AC = 10;
2685 Double_t minres[3]={20.,1,0.8};
2687 for(Int_t i=0;i<3;i++){
2689 fdriftA.ClearPoints();
2690 fdriftC.ClearPoints();
2691 fdriftAC.ClearPoints();
2700 for (Int_t itrack=0; itrack<338; ++itrack){
2701 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2702 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2704 //-- Exclude the tracks which has the biggest inclanation angle
2705 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2706 Int_t clustercounter=0;
2709 //-- exclude the low intensity tracks
2711 for (Int_t index=0; index<indexMax; ++index){
2713 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2714 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2715 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2717 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2719 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2724 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2726 if (itrack>335) indexMax=557568/2;
2727 for (Int_t index=0; index<indexMax; ++index){
2728 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2729 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2730 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2731 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2733 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2734 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2735 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2736 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2738 //cut if no track info available
2739 if (iltr->GetBundle()==0) continue;
2740 if (iR<133||mR<133) continue;
2741 if(mltr->fVecP2->GetMatrixArray()[index]>minres[i]) continue;
2743 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2744 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2746 //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2747 Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2749 if (iltr->GetSide()==0){
2750 fdriftA.AddPoint(xxx,mdrift,1);
2752 fdriftC.AddPoint(xxx,mdrift,1);
2754 // Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2755 Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2756 fdriftAC.AddPoint(xxx2,mdrift,1);
2759 }//end laser track loop
2769 fdriftA.GetParameters(fitA);
2770 fdriftC.GetParameters(fitC);
2771 fdriftAC.GetParameters(fitAC);
2773 //Parameters: 0 linear offset
2774 // 1 mean drift velocity correction factor
2775 // 2 relative global y gradient
2776 // 3 radial deformation
2777 // 4 IROC/OROC offset
2779 // FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2781 for (Int_t itrack=0; itrack<338; ++itrack){
2782 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2783 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2785 //-- Exclude the tracks which has the biggest inclanation angle
2786 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2787 Int_t clustercounter=0;
2790 //-- exclude the low intensity tracks
2792 for (Int_t index=0; index<indexMax; ++index){
2793 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2794 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2795 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2796 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2798 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2802 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2804 if (itrack>335) indexMax=557568/2;
2805 for (Int_t index=0; index<indexMax; ++index){
2806 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2807 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2808 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2809 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2811 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2812 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2813 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2814 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2816 //cut if no track info available
2817 if (iR<60||mR<60) continue;
2819 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2820 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2823 if (iltr->GetSide()==1) v=&fitC;
2824 // 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);
2825 Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2827 mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2836 //set statistics values
2838 npointsA= fdriftA.GetNpoints();
2839 if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2843 npointsC= fdriftC.GetNpoints();
2844 if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2848 npointsAC= fdriftAC.GetNpoints();
2849 if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2853 if (fStreamLevel>2){
2854 //laser track information
2855 (*GetDebugStreamer()) << "DriftV" <<
2857 "run=" << fRunNumber <<
2858 "time=" << timestamp <<
2859 "fitA.=" << &fitA <<
2860 "fitC.=" << &fitC <<
2861 "fitAC.=" << &fitAC <<
2871 //Parameters: 0 linear offset (global)
2872 // 1 mean drift velocity correction factor
2873 // 2 relative global y gradient
2874 // 3 radial deformation
2875 // 4 IROC/OROC offset
2876 // 5 linear offset relative A-C
2879 TGraphErrors *grA[7];
2880 TGraphErrors *grC[7];
2881 TGraphErrors *grAC[8];
2882 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_");
2883 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_");
2885 TObjArray *arrNames=names.Tokenize(";");
2886 TObjArray *arrNamesAC=namesAC.Tokenize(";");
2889 for (Int_t i=0; i<7; ++i){
2890 TString grName=arrNames->UncheckedAt(i)->GetName();
2892 grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2894 grA[i]=new TGraphErrors;
2895 grA[i]->SetName(grName.Data());
2896 grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2897 fArrFitGraphs->Add(grA[i]);
2899 // Int_t ipoint=grA[i]->GetN();
2901 grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2902 grA[i]->SetPointError(ipoint,60,.0001);
2903 if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2907 for (Int_t i=0; i<7; ++i){
2908 TString grName=arrNames->UncheckedAt(i)->GetName();
2910 grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2912 grC[i]=new TGraphErrors;
2913 grC[i]->SetName(grName.Data());
2914 grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2915 fArrFitGraphs->Add(grC[i]);
2917 // Int_t ipoint=grC[i]->GetN();
2919 grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2920 grC[i]->SetPointError(ipoint,60,.0001);
2921 if (i<4) grC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2925 for (Int_t i=0; i<8; ++i){
2926 TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2928 grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2930 grAC[i]=new TGraphErrors;
2931 grAC[i]->SetName(grName.Data());
2932 grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2933 fArrFitGraphs->Add(grAC[i]);
2935 // Int_t ipoint=grAC[i]->GetN();
2937 grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2938 grAC[i]->SetPointError(ipoint,60,.0001);
2939 if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2942 if (fDebugLevel>10){
2943 printf("A side fit parameters:\n");
2945 printf("\nC side fit parameters:\n");
2947 printf("\nAC side fit parameters:\n");
2954 //_____________________________________________________________________
2955 Double_t AliTPCCalibCE::SetBurstHnDrift()
2958 // Create a new THnSparse for the current burst
2959 // return the time of the current burst
2962 for(i=0; i<fTimeBursts.GetNrows(); ++i){
2963 if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2964 if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2965 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2966 return fTimeBursts(i);
2971 fArrHnDrift.AddAt(fHnDrift,i);
2972 fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
2977 //_____________________________________________________________________
2978 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
2981 // Write class to file
2982 // option can be specified in the dir option:
2984 // name=<objname>: the name of the calibration object in file will be <objname>
2985 // type=<type>: the saving type:
2986 // 0 - write the complte object
2987 // 1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
2988 // 2 - like 2, but in addition delete objects that will most probably not be used for calibration
2989 // 3 - store only calibration output, don't store the reference histograms
2990 // and THnSparse (requires Analyse called before)
2992 // NOTE: to read the object back, the ReadFromFile function should be used
2996 TString objName=GetName();
3000 TObjArray *arr=sDir.Tokenize(",");
3003 while ( (s=(TObjString*)next()) ){
3004 TString optString=s->GetString();
3005 optString.Remove(TString::kBoth,' ');
3006 if (optString.BeginsWith("name=")){
3007 objName=optString.ReplaceAll("name=","");
3009 if (optString.BeginsWith("type=")){
3010 optString.ReplaceAll("type=","");
3011 type=optString.Atoi();
3015 if (type==1||type==2) {
3016 //delete most probably not needed stuff
3017 //This requires Analyse to be called after reading the object from file
3018 fCalRocArrayT0.Delete();
3019 fCalRocArrayT0Err.Delete();
3020 fCalRocArrayQ.Delete();
3021 fCalRocArrayRMS.Delete();
3022 fCalRocArrayOutliers.Delete();
3024 if (type==2||type==3){
3025 fParamArrayEventPol1.Delete();
3026 fParamArrayEventPol2.Delete();
3029 TObjArray histoQArray(72);
3030 TObjArray histoT0Array(72);
3031 TObjArray histoRMSArray(72);
3032 TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3034 //save all large 2D histograms in separte pointers
3035 //to have a smaller memory print when saving the object
3036 if (type==1||type==2||type==3){
3037 for (Int_t i=0; i<72; ++i){
3038 histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3039 histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3040 histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3042 fHistoQArray.SetOwner(kFALSE);
3043 fHistoT0Array.SetOwner(kFALSE);
3044 fHistoRMSArray.SetOwner(kFALSE);
3045 fHistoQArray.Clear();
3046 fHistoT0Array.Clear();
3047 fHistoRMSArray.Clear();
3049 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3050 arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3052 fArrHnDrift.SetOwner(kFALSE);
3053 fArrHnDrift.Clear();
3057 TDirectory *backup = gDirectory;
3059 TFile f(filename,"recreate");
3060 this->Write(objName.Data());
3061 if (type==1||type==2) {
3062 histoQArray.Write("histoQArray",TObject::kSingleKey);
3063 histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3064 histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3065 arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3071 //move histograms back to the object
3072 if (type==1||type==2){
3073 for (Int_t i=0; i<72; ++i){
3074 fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3075 fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3076 fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3078 fHistoQArray.SetOwner(kTRUE);
3079 fHistoT0Array.SetOwner(kTRUE);
3080 fHistoRMSArray.SetOwner(kTRUE);
3082 for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3083 fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3085 fArrHnDrift.SetOwner(kTRUE);
3088 if ( backup ) backup->cd();
3090 //_____________________________________________________________________
3091 AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3094 // Read object from file
3095 // Handle properly if the histogram arrays were stored separately
3096 // call Analyse to make sure to have the calibration relevant information in the object
3100 if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3101 TList *l=f.GetListOfKeys();
3106 AliTPCCalibCE *ce=0x0;
3107 TObjArray *histoQArray=0x0;
3108 TObjArray *histoT0Array=0x0;
3109 TObjArray *histoRMSArray=0x0;
3110 TObjArray *arrHnDrift=0x0;
3112 while ( (key=(TKey*)next()) ){
3114 if ( o->IsA()==AliTPCCalibCE::Class() ){
3115 ce=(AliTPCCalibCE*)o;
3116 } else if ( o->IsA()==TObjArray::Class() ){
3117 TString name=key->GetName();
3118 if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3119 if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3120 if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3121 if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3126 //move histograms back to the object
3129 for (Int_t i=0; i<72; ++i){
3130 hist=(TH1*)histoQArray->UncheckedAt(i);
3131 if (hist) hist->SetDirectory(0x0);
3132 ce->fHistoQArray.AddAt(hist,i);
3134 ce->fHistoQArray.SetOwner(kTRUE);
3138 for (Int_t i=0; i<72; ++i){
3139 hist=(TH1*)histoT0Array->UncheckedAt(i);
3140 if (hist) hist->SetDirectory(0x0);
3141 ce->fHistoT0Array.AddAt(hist,i);
3143 ce->fHistoT0Array.SetOwner(kTRUE);
3147 for (Int_t i=0; i<72; ++i){
3148 hist=(TH1*)histoRMSArray->UncheckedAt(i);
3149 if (hist) hist->SetDirectory(0x0);
3150 ce->fHistoRMSArray.AddAt(hist,i);
3152 ce->fHistoRMSArray.SetOwner(kTRUE);
3156 for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3157 THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3158 ce->fArrHnDrift.AddAt(hSparse,i);