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>
284 #include <TSpectrum.h>
288 #include "AliRawReader.h"
289 #include "AliRawReaderRoot.h"
290 #include "AliRawReaderDate.h"
291 #include "AliRawEventHeaderBase.h"
292 #include "AliTPCRawStream.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),
384 // AliTPCSignal default constructor
386 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
390 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
391 for (Int_t i=0;i<14;++i){
395 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
397 //_____________________________________________________________________
398 AliTPCCalibCE::AliTPCCalibCE(const AliTPCCalibCE &sig) :
399 AliTPCCalibRawBase(sig),
400 fNbinsT0(sig.fNbinsT0),
401 fXminT0(sig.fXminT0),
402 fXmaxT0(sig.fXmaxT0),
403 fNbinsQ(sig.fNbinsQ),
406 fNbinsRMS(sig.fNbinsRMS),
407 fXminRMS(sig.fXminRMS),
408 fXmaxRMS(sig.fXmaxRMS),
409 fPeakDetMinus(sig.fPeakDetMinus),
410 fPeakDetPlus(sig.fPeakDetPlus),
411 fPeakIntMinus(sig.fPeakIntMinus),
412 fPeakIntPlus(sig.fPeakIntPlus),
413 fNoiseThresholdMax(sig.fNoiseThresholdMax),
414 fNoiseThresholdSum(sig.fNoiseThresholdSum),
415 fIsZeroSuppressed(sig.fIsZeroSuppressed),
418 fParam(new AliTPCParam),
424 fCalRocArrayT0Err(72),
427 fCalRocArrayOutliers(72),
431 fMeanT0rms(sig.fMeanT0rms),
432 fMeanQrms(sig.fMeanQrms),
433 fMeanRMSrms(sig.fMeanRMSrms),
435 fParamArrayEventPol1(72),
436 fParamArrayEventPol2(72),
437 fTMeanArrayEvent(72),
438 fQMeanArrayEvent(72),
439 fVEventTime(sig.fVEventTime),
440 fVEventNumber(sig.fVEventNumber),
441 fVTime0SideA(sig.fVTime0SideA),
442 fVTime0SideC(sig.fVTime0SideC),
445 fPadTimesArrayEvent(72),
447 fPadRMSArrayEvent(72),
448 fPadPedestalArrayEvent(72),
458 fVTime0OffsetCounter(72),
461 fCurrentCETimeRef(0),
462 fProcessOld(sig.fProcessOld),
463 fProcessNew(sig.fProcessNew),
464 fAnalyseNew(sig.fAnalyseNew),
472 // AliTPCSignal copy constructor
474 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
476 for (Int_t iSec = 0; iSec < 72; ++iSec){
477 const AliTPCCalROC *calQ = (AliTPCCalROC*)sig.fCalRocArrayQ.UncheckedAt(iSec);
478 const AliTPCCalROC *calT0 = (AliTPCCalROC*)sig.fCalRocArrayT0.UncheckedAt(iSec);
479 const AliTPCCalROC *calRMS = (AliTPCCalROC*)sig.fCalRocArrayRMS.UncheckedAt(iSec);
480 const AliTPCCalROC *calOut = (AliTPCCalROC*)sig.fCalRocArrayOutliers.UncheckedAt(iSec);
482 const TH2S *hQ = (TH2S*)sig.fHistoQArray.UncheckedAt(iSec);
483 const TH2S *hT0 = (TH2S*)sig.fHistoT0Array.UncheckedAt(iSec);
484 const TH2S *hRMS = (TH2S*)sig.fHistoRMSArray.UncheckedAt(iSec);
486 if ( calQ != 0x0 ) fCalRocArrayQ.AddAt(new AliTPCCalROC(*calQ), iSec);
487 if ( calT0 != 0x0 ) fCalRocArrayT0.AddAt(new AliTPCCalROC(*calT0), iSec);
488 if ( calRMS != 0x0 ) fCalRocArrayRMS.AddAt(new AliTPCCalROC(*calRMS), iSec);
489 if ( calOut != 0x0 ) fCalRocArrayOutliers.AddAt(new AliTPCCalROC(*calOut), iSec);
492 TH2S *hNew = new TH2S(*hQ);
493 hNew->SetDirectory(0);
494 fHistoQArray.AddAt(hNew,iSec);
497 TH2S *hNew = new TH2S(*hT0);
498 hNew->SetDirectory(0);
499 fHistoT0Array.AddAt(hNew,iSec);
502 TH2S *hNew = new TH2S(*hRMS);
503 hNew->SetDirectory(0);
504 fHistoRMSArray.AddAt(hNew,iSec);
508 //copy fit parameters event by event
510 for (Int_t iSec=0; iSec<72; ++iSec){
511 arr = (TObjArray*)sig.fParamArrayEventPol1.UncheckedAt(iSec);
513 TObjArray *arrEvents = new TObjArray(arr->GetSize());
514 fParamArrayEventPol1.AddAt(arrEvents, iSec);
515 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
516 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
517 arrEvents->AddAt(new TVectorD(*vec),iEvent);
520 arr = (TObjArray*)sig.fParamArrayEventPol2.UncheckedAt(iSec);
522 TObjArray *arrEvents = new TObjArray(arr->GetSize());
523 fParamArrayEventPol2.AddAt(arrEvents, iSec);
524 for (Int_t iEvent=0; iEvent<arr->GetSize(); ++iEvent)
525 if ( TVectorD *vec=(TVectorD*)arr->UncheckedAt(iEvent) )
526 arrEvents->AddAt(new TVectorD(*vec),iEvent);
529 TVectorF *vMeanTime = (TVectorF*)sig.fTMeanArrayEvent.UncheckedAt(iSec);
530 TVectorF *vMeanQ = (TVectorF*)sig.fQMeanArrayEvent.UncheckedAt(iSec);
532 fTMeanArrayEvent.AddAt(new TVectorF(*vMeanTime), iSec);
534 fQMeanArrayEvent.AddAt(new TVectorF(*vMeanQ), iSec);
538 fVEventTime.ResizeTo(sig.fVEventTime);
539 fVEventNumber.ResizeTo(sig.fVEventNumber);
540 fVEventTime.SetElements(sig.fVEventTime.GetMatrixArray());
541 fVEventNumber.SetElements(sig.fVEventNumber.GetMatrixArray());
545 for (Int_t i=0; i<sig.fArrHnDrift.GetEntries();++i){
546 TObject *o=sig.fArrHnDrift.UncheckedAt(i);
548 TObject *newo=o->Clone("fHnDrift");
549 fArrHnDrift.AddAt(newo,i);
550 if (sig.fHnDrift && o==sig.fHnDrift) fHnDrift=(THnSparseI*)newo;
554 for (Int_t i=0;i<sig.fTimeBursts.GetNrows();++i){
555 fTimeBursts[i]=sig.fTimeBursts[i];
558 for (Int_t i=0;i<14;++i){
559 fPeaks[i]=sig.fPeaks[i];
560 fPeakWidths[i]=sig.fPeakWidths[i];
562 if (sig.fArrFitGraphs) {
563 fArrFitGraphs=(TObjArray*)sig.fArrFitGraphs->Clone();
564 fArrFitGraphs->SetOwner();
567 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=sig.fBinsLastAna[i];
570 //_____________________________________________________________________
571 AliTPCCalibCE::AliTPCCalibCE(const TMap *config) :
572 AliTPCCalibRawBase(),
586 fNoiseThresholdMax(5.),
587 fNoiseThresholdSum(8.),
588 fIsZeroSuppressed(kFALSE),
591 fParam(new AliTPCParam),
597 fCalRocArrayT0Err(72),
600 fCalRocArrayOutliers(72),
608 fParamArrayEventPol1(72),
609 fParamArrayEventPol2(72),
610 fTMeanArrayEvent(72),
611 fQMeanArrayEvent(72),
618 fPadTimesArrayEvent(72),
620 fPadRMSArrayEvent(72),
621 fPadPedestalArrayEvent(72),
631 fVTime0OffsetCounter(72),
634 fCurrentCETimeRef(0),
645 // constructor which uses a tmap as input to set some specific parameters
647 SetNameTitle("AliTPCCalibCE","AliTPCCalibCE");
650 if (config->GetValue("FirstTimeBin")) fFirstTimeBin = ((TObjString*)config->GetValue("FirstTimeBin"))->GetString().Atoi();
651 if (config->GetValue("LastTimeBin")) fLastTimeBin = ((TObjString*)config->GetValue("LastTimeBin"))->GetString().Atoi();
652 if (config->GetValue("NbinsT0")) fNbinsT0 = ((TObjString*)config->GetValue("NbinsT0"))->GetString().Atoi();
653 if (config->GetValue("XminT0")) fXminT0 = ((TObjString*)config->GetValue("XminT0"))->GetString().Atof();
654 if (config->GetValue("XmaxT0")) fXmaxT0 = ((TObjString*)config->GetValue("XmaxT0"))->GetString().Atof();
655 if (config->GetValue("NbinsQ")) fNbinsQ = ((TObjString*)config->GetValue("NbinsQ"))->GetString().Atoi();
656 if (config->GetValue("XminQ")) fXminQ = ((TObjString*)config->GetValue("XminQ"))->GetString().Atof();
657 if (config->GetValue("XmaxQ")) fXmaxQ = ((TObjString*)config->GetValue("XmaxQ"))->GetString().Atof();
658 if (config->GetValue("NbinsRMS")) fNbinsRMS = ((TObjString*)config->GetValue("NbinsRMS"))->GetString().Atoi();
659 if (config->GetValue("XminRMS")) fXminRMS = ((TObjString*)config->GetValue("XminRMS"))->GetString().Atof();
660 if (config->GetValue("XmaxRMS")) fXmaxRMS = ((TObjString*)config->GetValue("XmaxRMS"))->GetString().Atof();
661 if (config->GetValue("PeakDetMinus")) fPeakDetMinus = ((TObjString*)config->GetValue("PeakDetMinus"))->GetString().Atoi();
662 if (config->GetValue("PeakDetPlus")) fPeakDetPlus = ((TObjString*)config->GetValue("PeakDetPlus"))->GetString().Atoi();
663 if (config->GetValue("PeakIntMinus")) fPeakIntMinus = ((TObjString*)config->GetValue("PeakIntMinus"))->GetString().Atoi();
664 if (config->GetValue("PeakIntPlus")) fPeakIntPlus = ((TObjString*)config->GetValue("PeakIntPlus"))->GetString().Atoi();
665 if (config->GetValue("NoiseThresholdMax")) fNoiseThresholdMax = ((TObjString*)config->GetValue("NoiseThresholdMax"))->GetString().Atof();
666 if (config->GetValue("NoiseThresholdSum")) fNoiseThresholdSum = ((TObjString*)config->GetValue("NoiseThresholdSum"))->GetString().Atof();
667 if (config->GetValue("IsZeroSuppressed")) fIsZeroSuppressed = (Bool_t)((TObjString*)config->GetValue("IsZeroSuppressed"))->GetString().Atoi();
668 if (config->GetValue("UseL1Phase")) fUseL1Phase = (Bool_t)((TObjString*)config->GetValue("UseL1Phase"))->GetString().Atoi();
669 if (config->GetValue("SecRejectRatio")) fSecRejectRatio = ((TObjString*)config->GetValue("SecRejectRatio"))->GetString().Atof();
671 if (config->GetValue("ProcessOld")) fProcessOld = (Bool_t)((TObjString*)config->GetValue("ProcessOld"))->GetString().Atoi();
672 if (config->GetValue("ProcessNew")) fProcessNew = (Bool_t)((TObjString*)config->GetValue("ProcessNew"))->GetString().Atoi();
673 if (config->GetValue("AnalyseNew")) fAnalyseNew = (Bool_t)((TObjString*)config->GetValue("AnalyseNew"))->GetString().Atoi();
675 for (Int_t i=0;i<1024;++i) fPadSignal[i]=0;
676 for (Int_t i=0;i<14;++i){
682 for (Int_t i=0; i<100; ++i) fBinsLastAna[i]=0;
685 //_____________________________________________________________________
686 AliTPCCalibCE& AliTPCCalibCE::operator = (const AliTPCCalibCE &source)
689 // assignment operator
691 if (&source == this) return *this;
692 new (this) AliTPCCalibCE(source);
696 //_____________________________________________________________________
697 AliTPCCalibCE::~AliTPCCalibCE()
703 fCalRocArrayT0.Delete();
704 fCalRocArrayT0Err.Delete();
705 fCalRocArrayQ.Delete();
706 fCalRocArrayRMS.Delete();
707 fCalRocArrayOutliers.Delete();
709 fHistoQArray.Delete();
710 fHistoT0Array.Delete();
711 fHistoRMSArray.Delete();
713 fHistoTmean.Delete();
715 fParamArrayEventPol1.Delete();
716 fParamArrayEventPol2.Delete();
717 fTMeanArrayEvent.Delete();
718 fQMeanArrayEvent.Delete();
720 fPadTimesArrayEvent.Delete();
721 fPadQArrayEvent.Delete();
722 fPadRMSArrayEvent.Delete();
723 fPadPedestalArrayEvent.Delete();
725 fArrHnDrift.SetOwner();
726 fArrHnDrift.Delete();
729 fArrFitGraphs->SetOwner();
730 delete fArrFitGraphs;
733 //_____________________________________________________________________
734 Int_t AliTPCCalibCE::Update(const Int_t icsector,
737 const Int_t icTimeBin,
738 const Float_t csignal)
741 // Signal filling methode on the fly pedestal and Time offset correction if necessary.
742 // no extra analysis necessary. Assumes knowledge of the signal shape!
743 // assumes that it is looped over consecutive time bins of one pad
746 if (!fProcessOld) return 0;
749 if (icRow<0) return 0;
750 if (icPad<0) return 0;
751 if (icTimeBin<0) return 0;
752 if ( (icTimeBin>fLastTimeBin) || (icTimeBin<fFirstTimeBin) ) return 0;
754 Int_t iChannel = fROC->GetRowIndexes(icsector)[icRow]+icPad; // global pad position in sector
756 //init first pad and sector in this event
757 if ( fCurrentChannel == -1 ) {
759 fCurrentChannel = iChannel;
760 fCurrentSector = icsector;
764 //process last pad if we change to a new one
765 if ( iChannel != fCurrentChannel ){
767 fLastSector=fCurrentSector;
768 fCurrentChannel = iChannel;
769 fCurrentSector = icsector;
773 //fill signals for current pad
774 fPadSignal[icTimeBin]=csignal;
775 if ( csignal > fMaxPadSignal ){
776 fMaxPadSignal = csignal;
777 fMaxTimeBin = icTimeBin;
782 //_____________________________________________________________________
783 void AliTPCCalibCE::ProcessBunch(const Int_t sector, const Int_t row, const Int_t pad,
784 const Int_t length, const UInt_t startTimeBin, const UShort_t* signal)
787 // new filling method to fill the THnSparse histogram
790 //only in new processing mode
791 if (!fProcessNew) return;
792 //don't use the IROCs and inner part of OROCs
793 if (sector<36) return;
795 //only bunches with reasonable length
796 if (length<3||length>10) return;
798 UShort_t timeBin = (UShort_t)startTimeBin;
799 //skip first laser layer
800 if (timeBin<280) return;
802 Double_t timeBurst=SetBurstHnDrift();
804 Int_t cePeak=((sector/18)%2)*7+6;
805 //after 1 event setup peak ranges
806 if (fEventInBunch==1 && fPeaks[cePeak]==0) {
808 fHnDrift->GetAxis(4)->SetRangeUser(timeBurst-2*60,timeBurst+2*60);
811 fHnDrift->GetAxis(4)->SetRangeUser(fHnDrift->GetAxis(4)->GetXmin(),fHnDrift->GetAxis(4)->GetXmax());
815 // After the first event only fill every 5th bin in a row with the CE information
817 if (fEventInBunch==0||(fPeaks[cePeak]>100&&TMath::Abs((Short_t)fPeaks[cePeak]-(Short_t)timeBin)<(Short_t)fPeakWidths[cePeak])){
824 if (!IsPeakInRange(timeBin+length/2,sector)) return;
826 Double_t x[kHnBinsDV]={(Double_t)sector,(Double_t)row,
827 (Double_t)padFill,(Double_t)timeBin,timeBurst};
829 for (Int_t iTimeBin = 0; iTimeBin<length; iTimeBin++){
830 Float_t sig=(Float_t)signal[iTimeBin];
831 // if (fPeaks[6]>900&&timeBin>(fPeaks[6]-20)&&sig<20) continue;
832 // if (fPeaks[6]>900&&timeBin<(fPeaks[6]-fPeakWidth[6])&&sig<5) continue;
834 fHnDrift->Fill(x,sig);
838 //_____________________________________________________________________
839 void AliTPCCalibCE::FindLaserLayers()
842 // Find the laser layer positoins
846 for (Int_t iside=0;iside<2;++iside){
848 //find CE signal position and width
849 fHnDrift->GetAxis(0)->SetRangeUser(36+iside*18,53+iside*18);
850 TH1D *hproj=fHnDrift->Projection(3);
851 hproj->GetXaxis()->SetRangeUser(700,1030);
852 Int_t maxbin=hproj->GetMaximumBin();
853 Double_t binc=hproj->GetBinCenter(maxbin);
854 hproj->GetXaxis()->SetRangeUser(binc-5,binc+5);
856 fPeaks[add+6]=(UShort_t)TMath::Nint(binc);
857 // fPeakWidths[4]=(UShort_t)TMath::Nint(4*hproj->GetRMS()+.5);
858 fPeakWidths[add+6]=7;
860 hproj->GetXaxis()->SetRangeUser(0,maxbin-10);
862 s.Search(hproj,2,"goff");
864 TMath::Sort(6,s.GetPositionX(),index,kFALSE);
865 for (Int_t i=0; i<6; ++i){
866 fPeaks[i+add]=(UShort_t)TMath::Nint(s.GetPositionX()[index[i]]);
867 fPeakWidths[i+add]=5;
872 // Int_t timepos=fPeaks[4]-2*fPeakWidths[4];
875 // for (Int_t i=3; i>=0; --i){
876 // hproj->GetXaxis()->SetRangeUser(timepos-width,timepos);
877 // fPeaks[i]=hproj->GetMaximumBin();
878 // fPeakWidths[i]=(UShort_t)TMath::Nint(10.);
880 // timepos=fPeaks[i]-width/2;
883 // for (Int_t i=add; i<add+7; ++i){
884 // printf("Peak: %u +- %u\n",fPeaks[i],fPeakWidths[i]);
886 //check width and reset peak if >100
887 // for (Int_t i=0; i<5; ++i){
888 // if (fPeakWidths[i]>100) {
898 //_____________________________________________________________________
899 void AliTPCCalibCE::FindPedestal(Float_t part)
902 // find pedestal and noise for the current pad. Use either database or
903 // truncated mean with part*100%
905 Bool_t noPedestal = kTRUE;
907 //use pedestal database if set
908 if (fPedestalTPC&&fPadNoiseTPC){
909 //only load new pedestals if the sector has changed
910 if ( fCurrentSector!=fLastSector ){
911 fPedestalROC = fPedestalTPC->GetCalROC(fCurrentSector);
912 fPadNoiseROC = fPadNoiseTPC->GetCalROC(fCurrentSector);
915 if ( fPedestalROC&&fPadNoiseROC ){
916 fPadPedestal = fPedestalROC->GetValue(fCurrentChannel)*(Float_t)(!fIsZeroSuppressed);
917 fPadNoise = fPadNoiseROC->GetValue(fCurrentChannel);
923 //if we are not running with pedestal database, or for the current sector there is no information
924 //available, calculate the pedestal and noise on the fly
928 if ( fIsZeroSuppressed ) return;
929 const Int_t kPedMax = 100; //maximum pedestal value
938 UShort_t histo[kPedMax];
939 memset(histo,0,kPedMax*sizeof(UShort_t));
941 //fill pedestal histogram
942 for (Int_t i=fFirstTimeBin; i<=fLastTimeBin; ++i){
943 padSignal = fPadSignal[i];
944 if (padSignal<=0) continue;
945 if (padSignal>max && i>10) {
949 if (padSignal>kPedMax-1) continue;
950 histo[int(padSignal+0.5)]++;
954 for (Int_t i=1; i<kPedMax; ++i){
955 if (count1<count0*0.5) median=i;
960 Float_t count=histo[median] ,mean=histo[median]*median, rms=histo[median]*median*median ;
962 for (Int_t idelta=1; idelta<10; ++idelta){
963 if (median-idelta<=0) continue;
964 if (median+idelta>kPedMax) continue;
965 if (count<part*count1){
966 count+=histo[median-idelta];
967 mean +=histo[median-idelta]*(median-idelta);
968 rms +=histo[median-idelta]*(median-idelta)*(median-idelta);
969 count+=histo[median+idelta];
970 mean +=histo[median+idelta]*(median+idelta);
971 rms +=histo[median+idelta]*(median+idelta)*(median+idelta);
976 rms = TMath::Sqrt(TMath::Abs(rms/count-mean*mean));
982 //_____________________________________________________________________
983 void AliTPCCalibCE::UpdateCETimeRef()
985 // Find the time reference of the last valid CE signal in sector
986 // for irocs of the A-Side the reference of the corresponging OROC is returned
987 // the reason are the non reflective bands on the A-Side, which make the reference very uncertain
988 if ( fLastSector == fCurrentSector ) return;
989 Int_t sector=fCurrentSector;
990 if ( sector < 18 ) sector+=36;
992 TVectorF *vtRef = GetTMeanEvents(sector);
993 if ( !vtRef ) return;
994 Int_t vtRefSize= vtRef->GetNrows();
995 if ( vtRefSize < fNevents+1 ) vtRef->ResizeTo(vtRefSize+100);
996 else vtRefSize=fNevents;
997 while ( (*vtRef)[vtRefSize]==0 && vtRefSize>=0 ) --vtRefSize;
998 fCurrentCETimeRef=(*vtRef)[vtRefSize];
999 AliDebug(3,Form("Sector: %02d - T0 ref: %.2f",fCurrentSector,fCurrentCETimeRef));
1001 //_____________________________________________________________________
1002 void AliTPCCalibCE::FindCESignal(TVectorD ¶m, Float_t &qSum, const TVectorF maxima)
1005 // Find position, signal width and height of the CE signal (last signal)
1006 // param[0] = Qmax, param[1] = mean time, param[2] = rms;
1007 // maxima: array of local maxima of the pad signal use the one closest to the mean CE position
1010 Float_t ceQmax =0, ceQsum=0, ceTime=0, ceRMS=0;
1012 Float_t ceSumThreshold = fNoiseThresholdSum*fPadNoise; // threshold for the signal sum
1013 const Int_t kCemin = fPeakIntMinus; // range for the analysis of the ce signal +- channels from the peak
1014 const Int_t kCemax = fPeakIntPlus;
1016 Float_t minDist = 25; //initial minimum distance betweek roc mean ce signal and pad ce signal
1018 // find maximum closest to the sector mean from the last event
1019 for ( Int_t imax=0; imax<maxima.GetNrows(); ++imax){
1020 // get sector mean of last event
1021 Float_t tmean = fCurrentCETimeRef;
1022 if ( TMath::Abs( tmean-maxima[imax] ) < minDist ) {
1023 minDist = tmean-maxima[imax];
1024 cemaxpos = (Int_t)maxima[imax];
1027 // printf("L1 phase TB: %f\n",GetL1PhaseTB());
1029 ceQmax = fPadSignal[cemaxpos]-fPadPedestal;
1030 for (Int_t i=cemaxpos-kCemin; i<=cemaxpos+kCemax; ++i){
1031 if ( (i>fFirstTimeBin) && (i<fLastTimeBin) ){
1032 Float_t signal = fPadSignal[i]-fPadPedestal;
1034 ceTime+=signal*(i+0.5);
1035 ceRMS +=signal*(i+0.5)*(i+0.5);
1041 if (ceQmax&&ceQsum>ceSumThreshold) {
1043 ceRMS = TMath::Sqrt(TMath::Abs(ceRMS/ceQsum-ceTime*ceTime));
1044 ceTime-=GetL1PhaseTB();
1045 fVTime0Offset.GetMatrixArray()[fCurrentSector]+=ceTime; // mean time for each sector
1046 fVTime0OffsetCounter.GetMatrixArray()[fCurrentSector]++;
1048 //Normalise Q to the 'cell-size': The wire density is the same in the IROC and OROC, therefore the
1049 // the pick-up signal should scale with the pad area. In addition
1050 // the signal should decrease with the wire distance (4mm in IROC, 6mm in OROC),
1051 // ratio 2/3. The pad area we express in cm2. We normalise the signal
1052 // to the OROC signal (factor 2/3 for the IROCs).
1053 Float_t norm = fParam->GetPadPitchWidth(fCurrentSector)*fParam->GetPadPitchLength(fCurrentSector,fCurrentRow);
1054 if ( fCurrentSector<fParam->GetNInnerSector() ) norm*=3./2.;
1057 fVMeanQ.GetMatrixArray()[fCurrentSector]+=ceQsum;
1058 fVMeanQCounter.GetMatrixArray()[fCurrentSector]++;
1070 //_____________________________________________________________________
1071 Bool_t AliTPCCalibCE::IsPeak(Int_t pos, Int_t tminus, Int_t tplus) const
1074 // Check if 'pos' is a Maximum. Consider 'tminus' timebins before
1075 // and 'tplus' timebins after 'pos'
1077 if ( (pos-tminus)<fFirstTimeBin || (pos+tplus)>fLastTimeBin ) return kFALSE;
1078 for (Int_t iTime = pos; iTime>pos-tminus; --iTime)
1079 if ( fPadSignal[iTime-1] >= fPadSignal[iTime] ) return kFALSE;
1080 for (Int_t iTime = pos, iTime2=pos; iTime<pos+tplus; ++iTime, ++iTime2){
1081 if ( (iTime==pos) && (fPadSignal[iTime+1]==fPadSignal[iTime]) ) // allow two timebins with same adc value
1083 if ( fPadSignal[iTime2+1] >= fPadSignal[iTime2] ) return kFALSE;
1087 //_____________________________________________________________________
1088 void AliTPCCalibCE::FindLocalMaxima(TVectorF &maxima)
1091 // Find local maxima on the pad signal and Histogram them
1093 Float_t ceThreshold = fNoiseThresholdMax*TMath::Max(fPadNoise,Float_t(1.)); // threshold for the signal
1096 for (Int_t i=fLastTimeBin-fPeakDetPlus+1; i>=fFirstTimeBin+fPeakDetMinus; --i){
1097 if ( (fPadSignal[i]-fPadPedestal)<ceThreshold ) continue;
1098 if (IsPeak(i,fPeakDetMinus,fPeakDetPlus) ){
1099 if (count<maxima.GetNrows()){
1100 maxima.GetMatrixArray()[count++]=i;
1101 GetHistoTmean(fCurrentSector,kTRUE)->Fill(i);
1102 i-=(fPeakDetMinus+fPeakDetPlus-1); // next peak cannot be at bin fPeakDetMinus+fPeakDetPlus-1
1107 //_____________________________________________________________________
1108 void AliTPCCalibCE::ProcessPad()
1111 // Process data of current pad
1115 TVectorF maxima(15); // the expected maximum number of maxima in the complete TPC should be 8 laser beam layers
1116 // + central electrode and possibly post peaks from the CE signal
1117 // however if we are on a high noise pad a lot more peaks due to the noise might occur
1118 FindLocalMaxima(maxima);
1119 if ( (fNevents == 0) || (fOldRunNumber!=fRunNumber) ) return; // return because we don't have Time0 info for the CE yet
1121 UpdateCETimeRef(); // update the time refenrence for the current sector
1122 if ( fCurrentCETimeRef<1e-30 ) return; //return if we don't have time 0 info, eg if only one side has laser
1125 FindCESignal(param, qSum, maxima);
1127 Double_t meanT = param[1];
1128 Double_t sigmaT = param[2];
1130 //Fill Event T0 counter
1131 (*GetPadTimesEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel] = meanT;
1134 GetHistoQ(fCurrentSector,kTRUE)->Fill( TMath::Sqrt(qSum), fCurrentChannel );
1136 //Fill RMS histogram
1137 GetHistoRMS(fCurrentSector,kTRUE)->Fill( sigmaT, fCurrentChannel );
1140 //Fill debugging info
1141 if ( GetStreamLevel()>0 ){
1142 (*GetPadPedestalEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=fPadPedestal;
1143 (*GetPadRMSEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=sigmaT;
1144 (*GetPadQEvent(fCurrentSector,kTRUE)).GetMatrixArray()[fCurrentChannel]=qSum;
1149 //_____________________________________________________________________
1150 void AliTPCCalibCE::EndEvent()
1152 // Process data of current pad
1153 // The Functions 'SetTimeStamp' and 'SetRunNumber' should be called
1154 // before the EndEvent function to set the event timestamp and number!!!
1155 // This is automatically done if the ProcessEvent(AliRawReader *rawReader)
1156 // function was called
1165 //check if last pad has allready been processed, if not do so
1166 if ( fMaxTimeBin>-1 ) ProcessPad();
1168 AliDebug(3, Form("EndEvent() - Start; Event: %05d", fNevents));
1171 TMatrixD dummy(3,3);
1172 // TVectorF vMeanTime(72);
1173 // TVectorF vMeanQ(72);
1174 AliTPCCalROC *calIroc=new AliTPCCalROC(0);
1175 AliTPCCalROC *calOroc=new AliTPCCalROC(36);
1177 //find mean time0 offset for side A and C
1178 //use only orocs due to the better statistics
1179 Double_t time0Side[2]; //time0 for side A:0 and C:1
1180 Double_t time0SideCount[2]; //time0 counter for side A:0 and C:1
1181 time0Side[0]=0;time0Side[1]=0;time0SideCount[0]=0;time0SideCount[1]=0;
1182 for ( Int_t iSec = 36; iSec<72; ++iSec ){
1183 time0Side[(iSec/18)%2] += fVTime0Offset.GetMatrixArray()[iSec];
1184 time0SideCount[(iSec/18)%2] += fVTime0OffsetCounter.GetMatrixArray()[iSec];
1186 if ( time0SideCount[0] >0 )
1187 time0Side[0]/=time0SideCount[0];
1188 if ( time0SideCount[1] >0 )
1189 time0Side[1]/=time0SideCount[1];
1190 // end find time0 offset
1191 AliDebug(3,Form("time0Side/time0SideCount: A=%.2f/%.2f, C=%.2f/%.2f",time0Side[0],time0SideCount[0],time0Side[1],time0SideCount[1]));
1193 //loop over all ROCs, fill CE Time histogram corrected for the mean Time0 of each ROC
1194 for ( Int_t iSec = 0; iSec<72; ++iSec ){
1195 AliDebug(4,Form("Processing sector '%02d'\n",iSec));
1196 //find median and then calculate the mean around it
1197 TH1S *hMeanT = GetHistoTmean(iSec); //histogram with local maxima position information
1198 if ( !hMeanT ) continue;
1199 //continue if not enough data is filled in the meanT histogram. This is the case if we do not have a laser event.
1200 if ( hMeanT->GetEffectiveEntries() < fROC->GetNChannels(iSec)*fSecRejectRatio ){
1202 AliDebug(3,Form("Skipping sec. '%02d': Not enough statistics\n",iSec));
1206 Double_t entries = hMeanT->GetEffectiveEntries();
1208 Short_t *arr = hMeanT->GetArray()+1;
1210 for ( ibin=0; ibin<hMeanT->GetNbinsX(); ++ibin){
1212 if ( sum>=(entries/2.) ) break;
1215 Int_t firstBin = fFirstTimeBin+ibin-delta;
1216 Int_t lastBin = fFirstTimeBin+ibin+delta;
1217 if ( firstBin<fFirstTimeBin ) firstBin=fFirstTimeBin;
1218 if ( lastBin>fLastTimeBin ) lastBin =fLastTimeBin;
1219 Float_t median =AliMathBase::GetCOG(arr+ibin-delta,2*delta,firstBin,lastBin);
1221 // check boundaries for ebye info of mean time
1222 TVectorF *vMeanTime=GetTMeanEvents(iSec,kTRUE);
1223 Int_t vSize=vMeanTime->GetNrows();
1224 if ( vSize < fNevents+1 ){
1225 vMeanTime->ResizeTo(vSize+100);
1228 // store mean time for the readout sides
1229 vSize=fVTime0SideA.GetNrows();
1230 if ( vSize < fNevents+1 ){
1231 fVTime0SideA.ResizeTo(vSize+100);
1232 fVTime0SideC.ResizeTo(vSize+100);
1234 fVTime0SideA.GetMatrixArray()[fNevents]=time0Side[0];
1235 fVTime0SideC.GetMatrixArray()[fNevents]=time0Side[1];
1237 vMeanTime->GetMatrixArray()[fNevents]=median;
1241 TVectorF *vTimes = GetPadTimesEvent(iSec);
1242 if ( !vTimes ) continue; //continue if no time information for this sector is available
1244 AliTPCCalROC calIrocOutliers(0);
1245 AliTPCCalROC calOrocOutliers(36);
1247 // calculate mean Q of the sector
1248 TVectorF *vMeanQ=GetQMeanEvents(iSec,kTRUE);
1249 vSize=vMeanQ->GetNrows();
1250 if ( vSize < fNevents+1 ){
1251 vMeanQ->ResizeTo(vSize+100);
1254 if ( fVMeanQCounter.GetMatrixArray()[iSec]>0 ) meanQ=fVMeanQ.GetMatrixArray()[iSec]/fVMeanQCounter.GetMatrixArray()[iSec];
1255 vMeanQ->GetMatrixArray()[fNevents]=meanQ;
1257 for ( UInt_t iChannel=0; iChannel<fROC->GetNChannels(iSec); ++iChannel ){
1258 Float_t time = (*vTimes).GetMatrixArray()[iChannel];
1260 //set values for temporary roc calibration class
1262 calIroc->SetValue(iChannel, time);
1263 if ( TMath::Abs(time) < 1e-30 ) calIrocOutliers.SetValue(iChannel,1);
1266 calOroc->SetValue(iChannel, time);
1267 if ( TMath::Abs(time) < 1e-30 ) calOrocOutliers.SetValue(iChannel,1);
1270 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) )
1271 // test that we really found the CE signal reliably
1272 if ( TMath::Abs(fVTime0SideA.GetMatrixArray()[fNevents-1]-time0Side[0])<.05)
1273 GetHistoT0(iSec,kTRUE)->Fill( time-time0Side[(iSec/18)%2],iChannel );
1277 //------------------------------- Debug start ------------------------------
1278 if ( GetStreamLevel()>0 ){
1279 TTreeSRedirector *streamer=GetDebugStreamer();
1285 Float_t q = (*GetPadQEvent(iSec))[iChannel];
1286 Float_t rms = (*GetPadRMSEvent(iSec))[iChannel];
1288 UInt_t channel=iChannel;
1291 while ( channel > (fROC->GetRowIndexes(sector)[row]+fROC->GetNPads(sector,row)-1) ) row++;
1292 pad = channel-fROC->GetRowIndexes(sector)[row];
1293 padc = pad-(fROC->GetNPads(sector,row)/2);
1295 // TH1F *h1 = new TH1F(Form("hSignalD%d.%d.%d",sector,row,pad),
1296 // Form("hSignalD%d.%d.%d",sector,row,pad),
1297 // fLastTimeBin-fFirstTimeBin,
1298 // fFirstTimeBin,fLastTimeBin);
1299 // h1->SetDirectory(0);
1301 // for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1302 // h1->Fill(i,fPadSignal(i));
1305 if (fVTime0OffsetCounter.GetMatrixArray()[iSec]>0)
1306 t0Sec = fVTime0Offset.GetMatrixArray()[iSec]/fVTime0OffsetCounter.GetMatrixArray()[iSec];
1307 Double_t t0Side = time0Side[(iSec/18)%2];
1308 (*streamer) << "DataPad" <<
1309 "Event=" << fNevents <<
1310 "RunNumber=" << fRunNumber <<
1311 "TimeStamp=" << fTimeStamp <<
1312 "Sector="<< sector <<
1316 "PadSec="<< channel <<
1317 "Time0Sec=" << t0Sec <<
1318 "Time0Side=" << t0Side <<
1322 "MeanQ=" << meanQ <<
1323 // "hist.=" << h1 <<
1329 //----------------------------- Debug end ------------------------------
1330 }// end channel loop
1333 //do fitting now only in debug mode
1334 if (GetDebugLevel()>0){
1335 TVectorD paramPol1(3);
1336 TVectorD paramPol2(6);
1337 TMatrixD matPol1(3,3);
1338 TMatrixD matPol2(6,6);
1342 if ( (fNevents>0) && (fOldRunNumber==fRunNumber) ){
1344 calIroc->GlobalFit(&calIrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1345 calIroc->GlobalFit(&calIrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1347 calOroc->GlobalFit(&calOrocOutliers,0,paramPol1,matPol1,chi2Pol1,0);
1348 calOroc->GlobalFit(&calOrocOutliers,0,paramPol2,matPol2,chi2Pol2,1);
1351 GetParamArrayPol1(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol1), fNevents);
1352 GetParamArrayPol2(iSec,kTRUE)->AddAtAndExpand(new TVectorD(paramPol2), fNevents);
1355 //------------------------------- Debug start ------------------------------
1356 if ( GetStreamLevel()>0 ){
1357 TTreeSRedirector *streamer=GetDebugStreamer();
1359 (*streamer) << "DataRoc" <<
1360 // "Event=" << fEvent <<
1361 "RunNumber=" << fRunNumber <<
1362 "TimeStamp=" << fTimeStamp <<
1364 "hMeanT.=" << hMeanT <<
1365 "median=" << median <<
1366 "paramPol1.=" << ¶mPol1 <<
1367 "paramPol2.=" << ¶mPol2 <<
1368 "matPol1.=" << &matPol1 <<
1369 "matPol2.=" << &matPol2 <<
1370 "chi2Pol1=" << chi2Pol1 <<
1371 "chi2Pol2=" << chi2Pol2 <<
1376 //------------------------------- Debug end ------------------------------
1379 //return if no sector has a valid mean time
1380 if ( nSecMeanT == 0 ) return;
1383 // fTMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanTime),fNevents);
1384 // fQMeanArrayEvent.AddAtAndExpand(new TVectorF(vMeanQ),fNevents);
1385 if ( fVEventTime.GetNrows() < fNevents+1 ) {
1386 fVEventTime.ResizeTo((Int_t)(fVEventTime.GetNrows()+100));
1387 fVEventNumber.ResizeTo((Int_t)(fVEventNumber.GetNrows()+100));
1389 fVEventTime.GetMatrixArray()[fNevents] = fTimeStamp;
1390 fVEventNumber.GetMatrixArray()[fNevents] = fEventId;
1393 if (fProcessNew) ++fEventInBunch;
1394 fOldRunNumber = fRunNumber;
1398 AliDebug(3, Form("EndEvent() - End; Event: %05d", fNevents));
1400 //_____________________________________________________________________
1401 TH2S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1402 Int_t nbinsY, Float_t ymin, Float_t ymax,
1403 const Char_t *type, Bool_t force)
1406 // return pointer to TH2S histogram of 'type'
1407 // if force is true create a new histogram if it doesn't exist allready
1409 if ( !force || arr->UncheckedAt(sector) )
1410 return (TH2S*)arr->UncheckedAt(sector);
1412 // if we are forced and histogram doesn't exist yet create it
1413 // new histogram with Q calib information. One value for each pad!
1414 TH2S* hist = new TH2S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1416 fROC->GetNChannels(sector),0,fROC->GetNChannels(sector));
1417 hist->SetDirectory(0);
1418 arr->AddAt(hist,sector);
1421 //_____________________________________________________________________
1422 TH2S* AliTPCCalibCE::GetHistoT0(Int_t sector, Bool_t force)
1425 // return pointer to T0 histogram
1426 // if force is true create a new histogram if it doesn't exist allready
1428 TObjArray *arr = &fHistoT0Array;
1429 return GetHisto(sector, arr, fNbinsT0, fXminT0, fXmaxT0, "T0", force);
1431 //_____________________________________________________________________
1432 TH2S* AliTPCCalibCE::GetHistoQ(Int_t sector, Bool_t force)
1435 // return pointer to Q histogram
1436 // if force is true create a new histogram if it doesn't exist allready
1438 TObjArray *arr = &fHistoQArray;
1439 return GetHisto(sector, arr, fNbinsQ, fXminQ, fXmaxQ, "Q", force);
1441 //_____________________________________________________________________
1442 TH2S* AliTPCCalibCE::GetHistoRMS(Int_t sector, Bool_t force)
1445 // return pointer to Q histogram
1446 // if force is true create a new histogram if it doesn't exist allready
1448 TObjArray *arr = &fHistoRMSArray;
1449 return GetHisto(sector, arr, fNbinsRMS, fXminRMS, fXmaxRMS, "RMS", force);
1451 //_____________________________________________________________________
1452 TH1S* AliTPCCalibCE::GetHisto(Int_t sector, TObjArray *arr,
1453 const Char_t *type, Bool_t force)
1456 // return pointer to TH1S histogram
1457 // if force is true create a new histogram if it doesn't exist allready
1459 if ( !force || arr->UncheckedAt(sector) )
1460 return (TH1S*)arr->UncheckedAt(sector);
1462 // if we are forced and histogram doesn't yes exist create it
1463 // new histogram with calib information. One value for each pad!
1464 TH1S* hist = new TH1S(Form("hCalib%s%.2d",type,sector),Form("%s calibration histogram sector %.2d",type,sector),
1465 fLastTimeBin-fFirstTimeBin,fFirstTimeBin,fLastTimeBin);
1466 hist->SetDirectory(0);
1467 arr->AddAt(hist,sector);
1470 //_____________________________________________________________________
1471 TH1S* AliTPCCalibCE::GetHistoTmean(Int_t sector, Bool_t force)
1474 // return pointer to Q histogram
1475 // if force is true create a new histogram if it doesn't exist allready
1477 TObjArray *arr = &fHistoTmean;
1478 return GetHisto(sector, arr, "LastTmean", force);
1480 //_____________________________________________________________________
1481 TVectorF* AliTPCCalibCE::GetVectSector(Int_t sector, TObjArray *arr, UInt_t size, Bool_t force) const
1484 // return pointer to Pad Info from 'arr' for the current event and sector
1485 // if force is true create it if it doesn't exist allready
1487 if ( !force || arr->UncheckedAt(sector) )
1488 return (TVectorF*)arr->UncheckedAt(sector);
1490 TVectorF *vect = new TVectorF(size);
1491 arr->AddAt(vect,sector);
1494 //_____________________________________________________________________
1495 TVectorF* AliTPCCalibCE::GetPadTimesEvent(Int_t sector, Bool_t force)
1498 // return pointer to Pad Times Array for the current event and sector
1499 // if force is true create it if it doesn't exist allready
1501 TObjArray *arr = &fPadTimesArrayEvent;
1502 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1504 //_____________________________________________________________________
1505 TVectorF* AliTPCCalibCE::GetPadQEvent(Int_t sector, Bool_t force)
1508 // return pointer to Pad Q Array for the current event and sector
1509 // if force is true create it if it doesn't exist allready
1510 // for debugging purposes only
1513 TObjArray *arr = &fPadQArrayEvent;
1514 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1516 //_____________________________________________________________________
1517 TVectorF* AliTPCCalibCE::GetPadRMSEvent(Int_t sector, Bool_t force)
1520 // return pointer to Pad RMS Array for the current event and sector
1521 // if force is true create it if it doesn't exist allready
1522 // for debugging purposes only
1524 TObjArray *arr = &fPadRMSArrayEvent;
1525 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1527 //_____________________________________________________________________
1528 TVectorF* AliTPCCalibCE::GetPadPedestalEvent(Int_t sector, Bool_t force)
1531 // return pointer to Pad RMS Array for the current event and sector
1532 // if force is true create it if it doesn't exist allready
1533 // for debugging purposes only
1535 TObjArray *arr = &fPadPedestalArrayEvent;
1536 return GetVectSector(sector,arr,fROC->GetNChannels(sector),force);
1538 //_____________________________________________________________________
1539 TVectorF* AliTPCCalibCE::GetTMeanEvents(Int_t sector, Bool_t force)
1542 // return pointer to the EbyE info of the mean arrival time for 'sector'
1543 // if force is true create it if it doesn't exist allready
1545 TObjArray *arr = &fTMeanArrayEvent;
1546 return GetVectSector(sector,arr,100,force);
1548 //_____________________________________________________________________
1549 TVectorF* AliTPCCalibCE::GetQMeanEvents(Int_t sector, Bool_t force)
1552 // return pointer to the EbyE info of the mean arrival time for 'sector'
1553 // if force is true create it if it doesn't exist allready
1555 TObjArray *arr = &fQMeanArrayEvent;
1556 return GetVectSector(sector,arr,100,force);
1558 //_____________________________________________________________________
1559 AliTPCCalROC* AliTPCCalibCE::GetCalRoc(Int_t sector, TObjArray* arr, Bool_t force) const
1562 // return pointer to ROC Calibration
1563 // if force is true create a new histogram if it doesn't exist allready
1565 if ( !force || arr->UncheckedAt(sector) )
1566 return (AliTPCCalROC*)arr->UncheckedAt(sector);
1568 // if we are forced and histogram doesn't yes exist create it
1570 // new AliTPCCalROC for T0 information. One value for each pad!
1571 AliTPCCalROC *croc = new AliTPCCalROC(sector);
1572 arr->AddAt(croc,sector);
1575 //_____________________________________________________________________
1576 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0(Int_t sector, Bool_t force)
1579 // return pointer to Time 0 ROC Calibration
1580 // if force is true create a new histogram if it doesn't exist allready
1582 TObjArray *arr = &fCalRocArrayT0;
1583 return GetCalRoc(sector, arr, force);
1585 //_____________________________________________________________________
1586 AliTPCCalROC* AliTPCCalibCE::GetCalRocT0Err(Int_t sector, Bool_t force)
1589 // return pointer to the error of Time 0 ROC Calibration
1590 // if force is true create a new histogram if it doesn't exist allready
1592 TObjArray *arr = &fCalRocArrayT0Err;
1593 return GetCalRoc(sector, arr, force);
1595 //_____________________________________________________________________
1596 AliTPCCalROC* AliTPCCalibCE::GetCalRocQ(Int_t sector, Bool_t force)
1599 // return pointer to T0 ROC Calibration
1600 // if force is true create a new histogram if it doesn't exist allready
1602 TObjArray *arr = &fCalRocArrayQ;
1603 return GetCalRoc(sector, arr, force);
1605 //_____________________________________________________________________
1606 AliTPCCalROC* AliTPCCalibCE::GetCalRocRMS(Int_t sector, Bool_t force)
1609 // return pointer to signal width ROC Calibration
1610 // if force is true create a new histogram if it doesn't exist allready
1612 TObjArray *arr = &fCalRocArrayRMS;
1613 return GetCalRoc(sector, arr, force);
1615 //_____________________________________________________________________
1616 AliTPCCalROC* AliTPCCalibCE::GetCalRocOutliers(Int_t sector, Bool_t force)
1619 // return pointer to Outliers
1620 // if force is true create a new histogram if it doesn't exist allready
1622 TObjArray *arr = &fCalRocArrayOutliers;
1623 return GetCalRoc(sector, arr, force);
1625 //_____________________________________________________________________
1626 TObjArray* AliTPCCalibCE::GetParamArray(Int_t sector, TObjArray* arr, Bool_t force) const
1629 // return pointer to TObjArray of fit parameters
1630 // if force is true create a new histogram if it doesn't exist allready
1632 if ( !force || arr->UncheckedAt(sector) )
1633 return (TObjArray*)arr->UncheckedAt(sector);
1635 // if we are forced and array doesn't yes exist create it
1637 // new TObjArray for parameters
1638 TObjArray *newArr = new TObjArray;
1639 arr->AddAt(newArr,sector);
1642 //_____________________________________________________________________
1643 TObjArray* AliTPCCalibCE::GetParamArrayPol1(Int_t sector, Bool_t force)
1646 // return pointer to TObjArray of fit parameters from plane fit
1647 // if force is true create a new histogram if it doesn't exist allready
1649 TObjArray *arr = &fParamArrayEventPol1;
1650 return GetParamArray(sector, arr, force);
1652 //_____________________________________________________________________
1653 TObjArray* AliTPCCalibCE::GetParamArrayPol2(Int_t sector, Bool_t force)
1656 // return pointer to TObjArray of fit parameters from parabola fit
1657 // if force is true create a new histogram if it doesn't exist allready
1659 TObjArray *arr = &fParamArrayEventPol2;
1660 return GetParamArray(sector, arr, force);
1663 //_____________________________________________________________________
1664 void AliTPCCalibCE::CreateDVhist()
1667 // Setup the THnSparse for the drift velocity determination
1671 //roc, row, pad, timebin, timestamp
1672 TTimeStamp begin(2010,01,01,0,0,0);
1673 TTimeStamp end(2035,01,01,0,0,0);
1674 Int_t nbinsTime=(end.GetSec()-begin.GetSec())/60; //Minutes resolution
1676 Int_t bins[kHnBinsDV] = { 72, 96, 140, 1030, nbinsTime};
1677 Double_t xmin[kHnBinsDV] = { 0., 0., 0., 0., (Double_t)begin.GetSec()};
1678 Double_t xmax[kHnBinsDV] = {72., 96., 140., 1030., (Double_t)end.GetSec()};
1680 fHnDrift=new THnSparseI("fHnDrift","Laser digits",kHnBinsDV, bins, xmin, xmax);
1681 fHnDrift->GetAxis(0)->SetNameTitle("ROC","Read-out chamber number");
1682 fHnDrift->GetAxis(1)->SetNameTitle("Row","Row number");
1683 fHnDrift->GetAxis(2)->SetNameTitle("Pad","Pad number");
1684 fHnDrift->GetAxis(3)->SetNameTitle("Timebin","Time bin [x100ns]");
1685 fHnDrift->GetAxis(4)->SetNameTitle("EventTime","Event time");
1689 //_____________________________________________________________________
1690 void AliTPCCalibCE::ResetEvent()
1693 // Reset global counters -- Should be called before each event is processed
1702 fPadTimesArrayEvent.Delete();
1703 fPadQArrayEvent.Delete();
1704 fPadRMSArrayEvent.Delete();
1705 fPadPedestalArrayEvent.Delete();
1707 for ( Int_t i=0; i<72; ++i ){
1708 fVTime0Offset.GetMatrixArray()[i]=0;
1709 fVTime0OffsetCounter.GetMatrixArray()[i]=0;
1710 fVMeanQ.GetMatrixArray()[i]=0;
1711 fVMeanQCounter.GetMatrixArray()[i]=0;
1714 //_____________________________________________________________________
1715 void AliTPCCalibCE::ResetPad()
1718 // Reset pad infos -- Should be called after a pad has been processed
1720 for (Int_t i=fFirstTimeBin; i<fLastTimeBin+1; ++i)
1727 //_____________________________________________________________________
1728 void AliTPCCalibCE::Merge(AliTPCCalibCE * const ce)
1731 // Merge ce to the current AliTPCCalibCE
1734 Int_t nCEevents = ce->GetNeventsProcessed();
1736 if (fProcessOld&&ce->fProcessOld){
1738 for (Int_t iSec=0; iSec<72; ++iSec){
1739 TH2S *hRefQmerge = ce->GetHistoQ(iSec);
1740 TH2S *hRefT0merge = ce->GetHistoT0(iSec);
1741 TH2S *hRefRMSmerge = ce->GetHistoRMS(iSec);
1745 TDirectory *dir = hRefQmerge->GetDirectory(); hRefQmerge->SetDirectory(0);
1746 TH2S *hRefQ = GetHistoQ(iSec);
1747 if ( hRefQ ) hRefQ->Add(hRefQmerge);
1749 TH2S *hist = new TH2S(*hRefQmerge);
1750 hist->SetDirectory(0);
1751 fHistoQArray.AddAt(hist, iSec);
1753 hRefQmerge->SetDirectory(dir);
1756 TDirectory *dir = hRefT0merge->GetDirectory(); hRefT0merge->SetDirectory(0);
1757 TH2S *hRefT0 = GetHistoT0(iSec);
1758 if ( hRefT0 ) hRefT0->Add(hRefT0merge);
1760 TH2S *hist = new TH2S(*hRefT0merge);
1761 hist->SetDirectory(0);
1762 fHistoT0Array.AddAt(hist, iSec);
1764 hRefT0merge->SetDirectory(dir);
1766 if ( hRefRMSmerge ){
1767 TDirectory *dir = hRefRMSmerge->GetDirectory(); hRefRMSmerge->SetDirectory(0);
1768 TH2S *hRefRMS = GetHistoRMS(iSec);
1769 if ( hRefRMS ) hRefRMS->Add(hRefRMSmerge);
1771 TH2S *hist = new TH2S(*hRefRMSmerge);
1772 hist->SetDirectory(0);
1773 fHistoRMSArray.AddAt(hist, iSec);
1775 hRefRMSmerge->SetDirectory(dir);
1780 // merge time information
1783 for (Int_t iSec=0; iSec<72; ++iSec){
1784 TObjArray *arrPol1CE = ce->GetParamArrayPol1(iSec);
1785 TObjArray *arrPol2CE = ce->GetParamArrayPol2(iSec);
1786 TVectorF *vMeanTimeCE = ce->GetTMeanEvents(iSec);
1787 TVectorF *vMeanQCE = ce->GetQMeanEvents(iSec);
1789 TObjArray *arrPol1 = 0x0;
1790 TObjArray *arrPol2 = 0x0;
1791 TVectorF *vMeanTime = 0x0;
1792 TVectorF *vMeanQ = 0x0;
1795 if ( arrPol1CE && arrPol2CE ){
1796 arrPol1 = GetParamArrayPol1(iSec,kTRUE);
1797 arrPol2 = GetParamArrayPol2(iSec,kTRUE);
1798 arrPol1->Expand(fNevents+nCEevents);
1799 arrPol2->Expand(fNevents+nCEevents);
1801 if ( vMeanTimeCE && vMeanQCE ){
1802 vMeanTime = GetTMeanEvents(iSec,kTRUE);
1803 vMeanQ = GetQMeanEvents(iSec,kTRUE);
1804 vMeanTime->ResizeTo(fNevents+nCEevents);
1805 vMeanQ->ResizeTo(fNevents+nCEevents);
1808 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1809 if ( arrPol1CE && arrPol2CE ){
1810 TVectorD *paramPol1 = (TVectorD*)(arrPol1CE->UncheckedAt(iEvent));
1811 TVectorD *paramPol2 = (TVectorD*)(arrPol2CE->UncheckedAt(iEvent));
1812 if ( paramPol1 && paramPol2 ){
1813 GetParamArrayPol1(iSec,kTRUE)->AddAt(new TVectorD(*paramPol1), fNevents+iEvent);
1814 GetParamArrayPol2(iSec,kTRUE)->AddAt(new TVectorD(*paramPol2), fNevents+iEvent);
1817 if ( vMeanTimeCE && vMeanQCE ){
1818 vMeanTime->GetMatrixArray()[fNevents+iEvent]=vMeanTimeCE->GetMatrixArray()[iEvent];
1819 vMeanQ->GetMatrixArray()[fNevents+iEvent]=vMeanQCE->GetMatrixArray()[iEvent];
1826 const TVectorD& eventTimes = ce->fVEventTime;
1827 const TVectorD& eventIds = ce->fVEventNumber;
1828 const TVectorF& time0SideA = ce->fVTime0SideA;
1829 const TVectorF& time0SideC = ce->fVTime0SideC;
1830 fVEventTime.ResizeTo(fNevents+nCEevents);
1831 fVEventNumber.ResizeTo(fNevents+nCEevents);
1832 fVTime0SideA.ResizeTo(fNevents+nCEevents);
1833 fVTime0SideC.ResizeTo(fNevents+nCEevents);
1835 for (Int_t iEvent=0; iEvent<nCEevents; ++iEvent){
1836 Double_t evTime = eventTimes.GetMatrixArray()[iEvent];
1837 Double_t evId = eventIds.GetMatrixArray()[iEvent];
1838 Float_t t0SideA = time0SideA.GetMatrixArray()[iEvent];
1839 Float_t t0SideC = time0SideC.GetMatrixArray()[iEvent];
1841 fVEventTime.GetMatrixArray()[fNevents+iEvent] = evTime;
1842 fVEventNumber.GetMatrixArray()[fNevents+iEvent] = evId;
1843 fVTime0SideA.GetMatrixArray()[fNevents+iEvent] = t0SideA;
1844 fVTime0SideC.GetMatrixArray()[fNevents+iEvent] = t0SideC;
1848 if (fProcessNew&&ce->fProcessNew) {
1849 if (fArrHnDrift.GetEntries() != ce->fArrHnDrift.GetEntries() ){
1850 AliError("Number of bursts in the instances to merge are different. No merging done!");
1852 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
1853 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
1854 THnSparseI *hce=(THnSparseI*)ce->fArrHnDrift.UncheckedAt(i);
1855 if (h && hce) h->Add(hce);
1856 else AliError(Form("AliTPCCalibCE::Merge - one THnSparse missing in burst %d",i));
1858 //TODO: What to do with fTimeBursts???
1862 fNevents+=nCEevents; //increase event counter
1865 //_____________________________________________________________________
1866 Long64_t AliTPCCalibCE::Merge(TCollection * const list)
1869 // Merge all objects of this type in list
1875 AliTPCCalibCE *ce=0;
1878 while ( (o=next()) ){
1879 ce=dynamic_cast<AliTPCCalibCE*>(o);
1889 //_____________________________________________________________________
1890 TGraph *AliTPCCalibCE::MakeGraphTimeCE(Int_t sector, Int_t xVariable, Int_t fitType, Int_t fitParameter)
1893 // Make graph from fit parameters of pol1 fit, pol2 fit, mean arrival time or mean Q for ROC 'sector'
1894 // or side (-1: A-Side, -2: C-Side)
1895 // xVariable: 0-event time, 1-event id, 2-internal event counter
1896 // fitType: 0-pol1 fit, 1-pol2 fit, 2-mean time, 3-mean Q
1897 // fitParameter: fit parameter ( 0-2 for pol1 ([0]+[1]*x+[2]*y),
1898 // 0-5 for pol2 ([0]+[1]*x+[2]*y+[3]*x*x+[4]*y*y+[5]*x*y),
1899 // not used for mean time and mean Q )
1900 // for an example see class description at the beginning
1903 TVectorD *xVar = 0x0;
1904 TObjArray *aType = 0x0;
1908 if ( (sector<-2) || (sector>71) ) return 0x0; //sector outside valid range
1909 if ( (xVariable<0) || (xVariable>2) ) return 0x0; //invalid x-variable
1910 if ( (fitType<0) || (fitType>3) ) return 0x0; //invalid fit type
1911 if ( sector>=0 && fitType==2 && !GetTMeanEvents(sector) ) return 0x0; //no mean time information available
1912 if ( sector>=0 && fitType==3 && !GetQMeanEvents(sector) ) return 0x0; //no mean charge information available
1913 if ( sector<0 && fitType!=2) return 0x0; //for side wise information only mean time is available
1917 if ( (fitParameter<0) || (fitParameter>2) ) return 0x0;
1918 aType = &fParamArrayEventPol1;
1919 if ( aType->At(sector)==0x0 ) return 0x0;
1921 else if ( fitType==1 ){
1922 if ( (fitParameter<0) || (fitParameter>5) ) return 0x0;
1923 aType = &fParamArrayEventPol2;
1924 if ( aType->At(sector)==0x0 ) return 0x0;
1928 if ( xVariable == 0 ) xVar = &fVEventTime;
1929 if ( xVariable == 1 ) xVar = &fVEventNumber;
1930 if ( xVariable == 2 ) {
1931 xVar = new TVectorD(fNevents);
1932 for ( Int_t i=0;i<fNevents; ++i) (*xVar)[i]=i;
1935 Double_t *x = new Double_t[fNevents];
1936 Double_t *y = new Double_t[fNevents];
1938 for (Int_t ievent =0; ievent<fNevents; ++ievent){
1940 TObjArray *events = (TObjArray*)(aType->At(sector));
1941 if ( events->GetSize()<=ievent ) break;
1942 TVectorD *v = (TVectorD*)(events->At(ievent));
1943 if ( (v!=0x0) && ((*xVar)[ievent]>0) ) { x[npoints]=(*xVar)[ievent]; y[npoints]=(*v)[fitParameter]; npoints++;}
1944 } else if (fitType == 2) {
1945 Double_t xValue=(*xVar)[ievent];
1947 if (sector>=0) yValue = (*GetTMeanEvents(sector))[ievent];
1948 else if (sector==-1) yValue=fVTime0SideA(ievent);
1949 else if (sector==-2) yValue=fVTime0SideC(ievent);
1950 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1951 }else if (fitType == 3) {
1952 Double_t xValue=(*xVar)[ievent];
1953 Double_t yValue=(*GetQMeanEvents(sector))[ievent];
1954 if ( yValue>0 && xValue>0 ) { x[npoints]=xValue; y[npoints]=yValue;npoints++;}
1958 TGraph *gr = new TGraph(npoints);
1959 //sort xVariable increasing
1960 Int_t *sortIndex = new Int_t[npoints];
1961 TMath::Sort(npoints,x,sortIndex, kFALSE);
1962 for (Int_t i=0;i<npoints;++i){
1963 gr->SetPoint(i,x[sortIndex[i]],y[sortIndex[i]]);
1967 if ( xVariable == 2 ) delete xVar;
1970 delete [] sortIndex;
1973 //_____________________________________________________________________
1974 void AliTPCCalibCE::Analyse()
1977 // Calculate calibration constants
1982 TVectorD paramT0(3);
1983 TVectorD paramRMS(3);
1984 TMatrixD dummy(3,3);
1986 Float_t channelCounter=0;
1991 for (Int_t iSec=0; iSec<72; ++iSec){
1992 TH2S *hT0 = GetHistoT0(iSec);
1993 if (!hT0 ) continue;
1995 AliTPCCalROC *rocQ = GetCalRocQ (iSec,kTRUE);
1996 AliTPCCalROC *rocT0 = GetCalRocT0 (iSec,kTRUE);
1997 AliTPCCalROC *rocT0Err = GetCalRocT0Err (iSec,kTRUE);
1998 AliTPCCalROC *rocRMS = GetCalRocRMS(iSec,kTRUE);
1999 AliTPCCalROC *rocOut = GetCalRocOutliers(iSec,kTRUE);
2001 TH2S *hQ = GetHistoQ(iSec);
2002 TH2S *hRMS = GetHistoRMS(iSec);
2004 Short_t *arrayhQ = hQ->GetArray();
2005 Short_t *arrayhT0 = hT0->GetArray();
2006 Short_t *arrayhRMS = hRMS->GetArray();
2008 UInt_t nChannels = fROC->GetNChannels(iSec);
2016 for (UInt_t iChannel=0; iChannel<nChannels; ++iChannel){
2019 Float_t cogTime0 = -1000;
2020 Float_t cogQ = -1000;
2021 Float_t cogRMS = -1000;
2027 Int_t offsetQ = (fNbinsQ+2)*(iChannel+1)+1;
2028 Int_t offsetT0 = (fNbinsT0+2)*(iChannel+1)+1;
2029 Int_t offsetRMS = (fNbinsRMS+2)*(iChannel+1)+1;
2031 cogQ = AliMathBase::GetCOG(arrayhQ+offsetQ,fNbinsQ,fXminQ,fXmaxQ,&rms);
2033 cogTime0 = AliMathBase::GetCOG(arrayhT0+offsetT0,fNbinsT0,fXminT0,fXmaxT0,&rmsT0);
2035 cogRMS = AliMathBase::GetCOG(arrayhRMS+offsetRMS,fNbinsRMS,fXminRMS,fXmaxRMS,&rms);
2040 //outlier specifications
2041 if ( (cogQ < ??) && (cogTime0 > ??) && (cogTime0<??) && ( cogRMS>??) ){
2048 rocQ->SetValue(iChannel, cogQ*cogQ);
2049 rocT0->SetValue(iChannel, cogTime0);
2050 rocT0Err->SetValue(iChannel, rmsT0);
2051 rocRMS->SetValue(iChannel, cogRMS);
2052 rocOut->SetValue(iChannel, cogOut);
2056 if ( GetStreamLevel() > 0 ){
2057 TTreeSRedirector *streamer=GetDebugStreamer();
2060 while ( iChannel > (fROC->GetRowIndexes(iSec)[row]+fROC->GetNPads(iSec,row)-1) ) row++;
2061 pad = iChannel-fROC->GetRowIndexes(iSec)[row];
2062 padc = pad-(fROC->GetNPads(iSec,row)/2);
2064 (*streamer) << "DataEnd" <<
2065 "Sector=" << iSec <<
2069 "PadSec=" << iChannel <<
2071 "T0=" << cogTime0 <<
2081 if ( channelCounter>0 ){
2082 fMeanT0rms/=channelCounter;
2083 fMeanQrms/=channelCounter;
2084 fMeanRMSrms/=channelCounter;
2086 // if ( fDebugStreamer ) fDebugStreamer->GetFile()->Write();
2087 // delete fDebugStreamer;
2088 // fDebugStreamer = 0x0;
2089 fVEventTime.ResizeTo(fNevents);
2090 fVEventNumber.ResizeTo(fNevents);
2091 fVTime0SideA.ResizeTo(fNevents);
2092 fVTime0SideC.ResizeTo(fNevents);
2095 if (fProcessNew&&fAnalyseNew){
2097 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries(); ++iburst){
2098 THnSparseI *h=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2099 h->GetAxis(4)->SetRangeUser(fTimeBursts[iburst]-60*5,fTimeBursts[iburst]+60*5);
2108 // New functions that also use the laser tracks
2113 //_____________________________________________________________________
2114 void AliTPCCalibCE::FindLocalMaxima(TObjArray * const arrObj, Double_t timestamp, Int_t burst)
2117 //Find the local maximums for the projections to each axis
2120 //find laser layer positoins
2121 fHnDrift->GetAxis(4)->SetRangeUser(timestamp-2*60,timestamp+2*60);
2123 THnSparse *hProj=fHnDrift;
2124 Double_t posCE[4]={0.,0.,0.,0.};
2125 Double_t widthCE[4]={0.,0.,0.,0.};
2127 // if(fPeaks[4]!=0){
2128 // find central electrode position once more, separately for IROC, OROC, A-, C-Side
2130 for (Int_t i=0; i<4; ++i){
2131 Int_t ce=(i/2>0)*7+6;
2132 hProj->GetAxis(0)->SetRangeUser(i*18,(i+1)*18-1);
2133 TH1 *h=fHnDrift->Projection(3);
2134 h->GetXaxis()->SetRangeUser(fPeaks[ce]-fPeakWidths[ce],fPeaks[ce]+fPeakWidths[ce]);
2135 Int_t nbinMax=h->GetMaximumBin();
2136 Double_t maximum=h->GetMaximum();
2137 // Double_t maxExpected=fNevents/fArrHnDrift->GetEntries()*556568./5./10.;
2138 // if (nbinMax<700||maximum<maxExpected) continue;
2139 Double_t xbinMax=h->GetBinCenter(nbinMax);
2140 TF1 fgaus("gaus","gaus",xbinMax-5,xbinMax+5);
2141 fgaus.SetParameters(maximum,xbinMax,2);
2142 fgaus.SetParLimits(1,xbinMax-5.,xbinMax+5.);
2143 fgaus.SetParLimits(2,0.2,4.);
2144 h->Fit(&fgaus,"RQN");
2145 // Double_t deltaX=4*fgaus.GetParameter(2);
2146 // xbinMax=fgaus.GetParameter(1);
2148 posCE[i]=fgaus.GetParameter(1);
2149 widthCE[i]=4*fgaus.GetParameter(2);
2150 hProj->GetAxis(0)->SetRangeUser(0,72);
2153 //Current drift velocity
2154 Float_t vdrift=2.61301900000000000e+06;//fParam->GetDriftV();
2155 // cout<<"vdrift="<<vdrift<<endl;
2157 AliDebug(5,Form("Timestamp %f - default drift velocity %f",timestamp,vdrift));
2158 //loop over all entries in the histogram
2160 for(Long64_t ichunk=0;ichunk<hProj->GetNbins();ichunk++){
2161 //get entry position and content
2162 Double_t adc=hProj->GetBinContent(ichunk,coord);
2165 Int_t sector = coord[0]-1;
2166 Int_t row = coord[1]-1;
2167 Int_t pad = coord[2]-1;
2168 Int_t timeBin= coord[3]-1;
2169 Double_t time = fHnDrift->GetAxis(4)->GetBinCenter(coord[4]);
2170 Int_t side = (sector/18)%2;
2172 // fPeaks[4]=(UInt_t)posCE[sector/18];
2173 // fPeakWidths[4]=(UInt_t)widthCE[sector/18];
2176 if (time<timestamp-120||time>timestamp+120) continue; //window of +- 2min
2177 if (adc < 5 ) continue;
2178 if (IsEdgePad(sector,row,pad)) continue;
2179 // if (!IsPeakInRange(timeBin)) continue;
2180 // if (isCE&&((row%2)||(row%2)||(sector%2))) continue;
2181 // if (isCE&&(sector!=0)) continue;
2183 Int_t padmin=-2, padmax=2;
2184 Int_t timemin=-2, timemax=2;
2185 Int_t minsumperpad=2;
2186 //CE or laser tracks
2188 if (TMath::Abs((Short_t)timeBin-(Short_t)posCE[sector/18])<(Short_t)widthCE[sector/18]) {
2197 // Find local maximum and cogs
2199 Bool_t isMaximum=kTRUE;
2200 Float_t totalmass=0, tbcm=0, padcm=0, rmstb=0, rmspad=0;
2201 Double_t cogY=0, rmsY=0;
2204 // for position calculation use
2205 for(Int_t ipad=padmin;ipad<=padmax;++ipad){
2207 fROC->GetPositionLocal(sector,row,pad+ipad,lxyz);
2209 for(Int_t itime=timemin;itime<=timemax;++itime){
2211 Int_t a[5]={coord[0],coord[1],coord[2]+ipad,coord[3]+itime,coord[4]};
2212 Double_t val=hProj->GetBinContent(a);
2215 tbcm +=(timeBin+itime)*val;
2216 padcm+=(pad+ipad)*val;
2219 rmstb +=(timeBin+itime)*(timeBin+itime)*val;
2220 rmspad+=(pad+ipad)*(pad+ipad)*val;
2221 rmsY +=lxyz[1]*lxyz[1]*val;
2229 if (!isMaximum) break;
2232 if (!isMaximum||!npart) continue;
2233 if (totalmass<npart*minsumperpad) continue;
2234 if (!isCE&&rmspad<.1) continue; //most probably noise, since signal only in one pad,
2235 //for CE we want only one pad by construction
2245 rmstb=TMath::Sqrt(TMath::Abs(tbcm*tbcm-rmstb));
2246 rmspad=TMath::Sqrt(TMath::Abs(padcm*padcm-rmspad));
2247 rmsY=TMath::Sqrt(TMath::Abs(cogY*cogY-rmsY));
2249 Int_t cog=TMath::Nint(padcm);
2251 // timebin --> z position
2252 Float_t zlength=fParam->GetZLength(side);
2253 // Float_t timePos=tbcm+(1000-fPeaks[4])
2254 // drift velocity is in m/s we would like to have cm/100ns, so 100cm/(10^7*100ns)
2255 Double_t gz=(zlength-(tbcm*vdrift*1.e-7))*TMath::Power(-1,side);
2257 // local to global transformation--> x and y positions
2259 fROC->GetPositionLocal(sector,row,pad,padlxyz);
2261 Double_t gxyz[3]={padlxyz[0],cogY,gz};
2262 Double_t lxyz[3]={padlxyz[0],cogY,gz};
2263 Double_t igxyz[3]={0,0,0};
2265 t1.RotatedGlobal2Global(sector,gxyz);
2271 //find track id and index of the position in the track (row)
2274 index=row+(sector>35)*fROC->GetNRows(0);
2275 trackID=FindLaserTrackID(sector,index,gxyz,mindist,lxyz,trackID2);
2277 trackID=336+((sector/18)%2);
2278 index= fROC->GetRowIndexes(sector)[row]+pad; // global pad position in sector
2280 index+=(sector%18)*fROC->GetNChannels(sector);
2282 index+=18*fROC->GetNChannels(0);
2283 index+=(sector%18)*fROC->GetNChannels(sector);
2285 //TODO: find out about the multiple peaks in the CE
2286 // mindist=TMath::Abs(fPeaks[4]-tbcm);
2290 // fill track vectors
2292 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arrObj->UncheckedAt(trackID);
2293 Double_t oldMinDist=ltr->fVecPhi->GetMatrixArray()[index];
2296 // travel time effect of light includes
2298 Double_t raylength=ltr->GetRayLength();
2299 Double_t globmir[3]={ltr->Xv(),ltr->Yv(),ltr->Zv()};
2300 ltr->GetXYZ(globmir);
2303 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2304 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2305 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2308 gxyz[2]=gxyz[2]-(TMath::Sqrt((gxyz[0]-globmir[0])*(gxyz[0]-globmir[0])
2309 +(gxyz[1]-globmir[1])*(gxyz[1]-globmir[1])
2310 +(gxyz[2]-globmir[2])*(gxyz[2]-globmir[2])+raylength))*vdrift*TMath::Power(10.,-6.)/30000;
2314 if (TMath::Abs(oldMinDist)<1.e-20||oldMinDist>mindist){
2315 ltr->fVecSec->GetMatrixArray()[index]=sector;
2316 ltr->fVecP2->GetMatrixArray()[index]=0;
2317 ltr->fVecPhi->GetMatrixArray()[index]=mindist;
2318 ltr->fVecGX->GetMatrixArray()[index]=gxyz[0];
2319 ltr->fVecGY->GetMatrixArray()[index]=gxyz[1];
2320 ltr->fVecGZ->GetMatrixArray()[index]=gxyz[2];
2321 ltr->fVecLX->GetMatrixArray()[index]=lxyz[0];
2322 ltr->fVecLY->GetMatrixArray()[index]=lxyz[1];
2323 ltr->fVecLZ->GetMatrixArray()[index]=lxyz[2];
2324 // ltr->SetUniqueID((UInt_t)(mindist*10000)); //distance in um
2326 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2327 ltr=(AliTPCLaserTrack*)arr->UncheckedAt(trackID);
2328 igxyz[0]=ltr->fVecGX->GetMatrixArray()[row];
2329 igxyz[1]=ltr->fVecGY->GetMatrixArray()[row];
2330 igxyz[2]=ltr->fVecGZ->GetMatrixArray()[row];
2334 if (fStreamLevel>4){
2335 (*GetDebugStreamer()) << "clusters" <<
2336 "run=" << fRunNumber <<
2337 "timestamp=" << timestamp <<
2338 "burst=" << burst <<
2344 "timebin=" << timeBin <<
2345 "cogCE=" << posCE[sector/18] <<
2346 "withCE=" << widthCE[sector/18] <<
2347 "index=" << index <<
2349 "padcm=" << padcm <<
2350 "rmspad=" << rmspad <<
2353 "rmstb=" << rmstb <<
2357 "lx=" << padlxyz[0]<<
2359 "lypad=" << padlxyz[1]<<
2366 "igx=" << igxyz[0] <<
2367 "igy=" << igxyz[1] <<
2368 "igz=" << igxyz[2] <<
2370 "mind=" << mindist <<
2372 "trackid=" << trackID <<
2373 "trackid2=" << trackID2 <<
2374 "npart=" << npart <<
2376 } // end stream levelmgz.fElements
2382 //_____________________________________________________________________
2383 void AliTPCCalibCE::AnalyseTrack()
2386 // Analyse the tracks
2390 AliTPCLaserTrack::LoadTracks();
2391 // AliTPCParam *param=0x0;
2393 // AliCDBManager *man=AliCDBManager::Instance();
2394 // if (man->GetDefaultStorage()){
2395 // AliCDBEntry *entry=man->Get("TPC/Calib/Parameters",fRunNumber);
2397 // entry->SetOwner(kTRUE);
2398 // param = (AliTPCParam*)(entry->GetObject()->Clone());
2402 // if (fParam) delete fParam;
2405 // AliError("Could not get updated AliTPCParam from OCDB!!!");
2408 //Measured and ideal laser tracks
2409 TObjArray* arrMeasured = SetupMeasured();
2410 TObjArray* arrIdeal = AliTPCLaserTrack::GetTracks();
2411 AddCEtoIdeal(arrIdeal);
2413 //find bursts and loop over them
2414 for (Int_t iburst=0; iburst<fArrHnDrift.GetEntries();++iburst){
2415 Double_t timestamp=fTimeBursts[iburst];
2416 AliDebug(5,Form("Burst: %d (%f)",iburst,timestamp));
2417 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(iburst);
2418 if (!fHnDrift) continue;
2419 UInt_t entries=(UInt_t)fHnDrift->GetEntries();
2420 if (fBinsLastAna[iburst]>=entries) continue; //already analysed!!!
2421 fBinsLastAna[iburst]=entries;
2423 for (Int_t iDim=0; iDim<fHnDrift->GetNdimensions(); ++iDim) fHnDrift->GetAxis(iDim)->SetRange(0,0);
2424 // if (iburst==0) FindLaserLayers();
2426 //reset laser tracks
2427 ResetMeasured(arrMeasured);
2429 //find clusters and associate to the tracks
2430 FindLocalMaxima(arrMeasured, timestamp, iburst);
2432 //calculate drift velocity
2433 CalculateDV(arrIdeal,arrMeasured,iburst);
2435 //Dump information to file if requested
2436 if (fStreamLevel>2){
2437 //printf("make tree\n");
2438 //laser track information
2440 for (Int_t itrack=0; itrack<338; ++itrack){
2441 TObject *iltr=arrIdeal->UncheckedAt(itrack);
2442 TObject *mltr=arrMeasured->UncheckedAt(itrack);
2443 (*GetDebugStreamer()) << "tracks" <<
2444 "run=" << fRunNumber <<
2445 "time=" << timestamp <<
2446 "burst="<< iburst <<
2453 if (fStreamLevel>0) GetDebugStreamer()->GetFile()->Write();
2456 //_____________________________________________________________________
2457 Int_t AliTPCCalibCE::FindLaserTrackID(Int_t sector,Int_t row, const Double_t *peakpos,Double_t &mindist,
2458 const Double_t *peakposloc, Int_t &itrackMin2)
2461 // Find the tracks, which are closest to the ideal tracks, from clusters closest to the ideal tracks
2465 TObjArray *arr=AliTPCLaserTrack::GetTracks();
2466 TVector3 vP(peakpos[0],peakpos[1],peakpos[2]);
2471 Int_t lastbeam=336/2;
2472 if ( (sector/18)%2 ) {
2479 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2480 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2481 // if (ltr->GetVecSec()->GetMatrixArray()[row]!=sector) continue;
2482 vSt.SetXYZ(ltr->GetX(),ltr->GetY(),ltr->GetZ());
2483 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2484 if (TMath::Abs(deltaZ)>40) continue;
2485 vDir.SetMagThetaPhi(1,ltr->Theta(),TMath::ASin(ltr->GetSnp()));
2486 vSt.RotateZ(ltr->GetAlpha());
2487 vDir.RotateZ(ltr->GetAlpha());
2489 Double_t dist=(vDir.Cross(vSt-vP)).Mag()/vDir.Mag();
2498 Float_t mindist2=10;
2499 for (Int_t itrack=firstbeam; itrack<lastbeam; ++itrack){
2500 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->At(itrack); //get the track
2501 if ((ltr->fVecSec->GetMatrixArray())[row]!=sector) continue;
2503 Double_t deltaZ=ltr->GetZ()-peakpos[2];
2504 if (TMath::Abs(deltaZ)>40) continue;
2506 Double_t dist=TMath::Abs((ltr->fVecLY->GetMatrixArray())[row]-peakposloc[1]);
2507 if (dist>1) continue;
2519 //_____________________________________________________________________
2520 Bool_t AliTPCCalibCE::IsEdgePad(Int_t sector, Int_t row, Int_t pad) const
2523 // return true if pad is on the edge of a row
2526 if ( pad == edge1 ) return kTRUE;
2527 Int_t edge2 = fROC->GetNPads(sector,row)-1;
2528 if ( pad == edge2 ) return kTRUE;
2533 //_____________________________________________________________________
2534 TObjArray* AliTPCCalibCE::SetupMeasured()
2537 // setup array of measured laser tracks and CE
2540 TObjArray *arrIdeal = AliTPCLaserTrack::GetTracks();
2541 TObjArray *arrMeasured = new TObjArray(338);
2542 arrMeasured->SetOwner();
2543 for(Int_t itrack=0;itrack<336;itrack++){
2544 arrMeasured->AddAt(new AliTPCLaserTrack(*((AliTPCLaserTrack*)arrIdeal->At(itrack))),itrack);
2548 AliTPCLaserTrack *ltrce=new AliTPCLaserTrack;
2551 ltrce->fVecSec=new TVectorD(557568/2);
2552 ltrce->fVecP2=new TVectorD(557568/2);
2553 ltrce->fVecPhi=new TVectorD(557568/2);
2554 ltrce->fVecGX=new TVectorD(557568/2);
2555 ltrce->fVecGY=new TVectorD(557568/2);
2556 ltrce->fVecGZ=new TVectorD(557568/2);
2557 ltrce->fVecLX=new TVectorD(557568/2);
2558 ltrce->fVecLY=new TVectorD(557568/2);
2559 ltrce->fVecLZ=new TVectorD(557568/2);
2561 arrMeasured->AddAt(ltrce,336); //CE A-Side
2563 ltrce=new AliTPCLaserTrack;
2566 ltrce->fVecSec=new TVectorD(557568/2);
2567 ltrce->fVecP2=new TVectorD(557568/2);
2568 ltrce->fVecPhi=new TVectorD(557568/2);
2569 ltrce->fVecGX=new TVectorD(557568/2);
2570 ltrce->fVecGY=new TVectorD(557568/2);
2571 ltrce->fVecGZ=new TVectorD(557568/2);
2572 ltrce->fVecLX=new TVectorD(557568/2);
2573 ltrce->fVecLY=new TVectorD(557568/2);
2574 ltrce->fVecLZ=new TVectorD(557568/2);
2575 arrMeasured->AddAt(ltrce,337); //CE C-Side
2580 //_____________________________________________________________________
2581 void AliTPCCalibCE::ResetMeasured(TObjArray * const arr)
2584 // reset array of measured laser tracks and CE
2586 for(Int_t itrack=0;itrack<338;itrack++){
2587 AliTPCLaserTrack *ltr=(AliTPCLaserTrack*)arr->UncheckedAt(itrack);
2588 ltr->fVecSec->Zero();
2589 ltr->fVecP2->Zero();
2590 ltr->fVecPhi->Zero();
2591 ltr->fVecGX->Zero();
2592 ltr->fVecGY->Zero();
2593 ltr->fVecGZ->Zero();
2594 ltr->fVecLX->Zero();
2595 ltr->fVecLY->Zero();
2596 ltr->fVecLZ->Zero();
2600 //_____________________________________________________________________
2601 void AliTPCCalibCE::AddCEtoIdeal(TObjArray *arr)
2604 // Add ideal CE positions to the ideal track data
2609 AliTPCLaserTrack *ltrceA=new AliTPCLaserTrack;
2612 ltrceA->fVecSec=new TVectorD(557568/2);
2613 ltrceA->fVecP2=new TVectorD(557568/2);
2614 ltrceA->fVecPhi=new TVectorD(557568/2);
2615 ltrceA->fVecGX=new TVectorD(557568/2);
2616 ltrceA->fVecGY=new TVectorD(557568/2);
2617 ltrceA->fVecGZ=new TVectorD(557568/2);
2618 ltrceA->fVecLX=new TVectorD(557568/2);
2619 ltrceA->fVecLY=new TVectorD(557568/2);
2620 ltrceA->fVecLZ=new TVectorD(557568/2);
2621 arr->AddAt(ltrceA,336); //CE A-Side
2623 AliTPCLaserTrack *ltrceC=new AliTPCLaserTrack;
2626 ltrceC->fVecSec=new TVectorD(557568/2);
2627 ltrceC->fVecP2=new TVectorD(557568/2);
2628 ltrceC->fVecPhi=new TVectorD(557568/2);
2629 ltrceC->fVecGX=new TVectorD(557568/2);
2630 ltrceC->fVecGY=new TVectorD(557568/2);
2631 ltrceC->fVecGZ=new TVectorD(557568/2);
2632 ltrceC->fVecLX=new TVectorD(557568/2);
2633 ltrceC->fVecLY=new TVectorD(557568/2);
2634 ltrceC->fVecLZ=new TVectorD(557568/2);
2635 arr->AddAt(ltrceC,337); //CE C-Side
2637 //Calculate ideal positoins
2640 Int_t channelSideA=0;
2641 Int_t channelSideC=0;
2642 Int_t channelSide=0;
2643 AliTPCLaserTrack *ltrce=0x0;
2645 for (Int_t isector=0; isector<72; ++isector){
2646 Int_t side=((isector/18)%2);
2647 for (UInt_t irow=0;irow<fROC->GetNRows(isector);++irow){
2648 for (UInt_t ipad=0;ipad<fROC->GetNPads(isector,irow);++ipad){
2649 fROC->GetPositionGlobal(isector,irow,ipad,gxyz);
2650 fROC->GetPositionLocal(isector,irow,ipad,lxyz);
2653 channelSide=channelSideA;
2656 channelSide=channelSideC;
2659 ltrce->fVecSec->GetMatrixArray()[channelSide]=isector;
2660 ltrce->fVecP2->GetMatrixArray()[channelSide]=0;
2661 ltrce->fVecPhi->GetMatrixArray()[channelSide]=0;
2662 ltrce->fVecGX->GetMatrixArray()[channelSide]=gxyz[0];
2663 ltrce->fVecGY->GetMatrixArray()[channelSide]=gxyz[1];
2664 // ltrce->fVecGZ->GetMatrixArray()[channelSide]=-1;
2665 ltrce->fVecLX->GetMatrixArray()[channelSide]=lxyz[0];
2666 ltrce->fVecLY->GetMatrixArray()[channelSide]=lxyz[1];
2667 // ltrce->fVecLZ->GetMatrixArray()[channelSide]=-1;
2670 ltrce->fVecGZ->GetMatrixArray()[channelSide]=-0.335;
2671 ltrce->fVecLZ->GetMatrixArray()[channelSide]=-0.335;
2675 ltrce->fVecGZ->GetMatrixArray()[channelSide]=0.15;
2676 ltrce->fVecLZ->GetMatrixArray()[channelSide]=0.15;
2686 //_____________________________________________________________________
2687 void AliTPCCalibCE::CalculateDV(TObjArray * const arrIdeal, TObjArray * const arrMeasured, Int_t burst)
2690 // calculate the drift velocity from the reconstructed clusters associated
2691 // to the ideal laser tracks
2692 // use two different fit scenarios: Separate fit for A- and C-Side
2693 // Common fit for A- and C-Side
2696 if (!fArrFitGraphs){
2697 fArrFitGraphs=new TObjArray;
2700 // static TLinearFitter fdriftA(5,"hyp4");
2701 // static TLinearFitter fdriftC(5,"hyp4");
2702 // static TLinearFitter fdriftAC(6,"hyp5");
2703 Double_t timestamp=fTimeBursts[burst];
2705 static TLinearFitter fdriftA(4,"hyp3");
2706 static TLinearFitter fdriftC(4,"hyp3");
2707 static TLinearFitter fdriftAC(5,"hyp4");
2708 TVectorD fitA(7),fitC(7),fitAC(8); //fit values+chi2+npoints
2712 Float_t chi2AC = 10;
2717 Double_t minres[3]={20.,1,0.8};
2719 for(Int_t i=0;i<3;i++){
2721 fdriftA.ClearPoints();
2722 fdriftC.ClearPoints();
2723 fdriftAC.ClearPoints();
2732 for (Int_t itrack=0; itrack<338; ++itrack){
2733 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2734 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2736 //-- Exclude the tracks which has the biggest inclanation angle
2737 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2738 Int_t clustercounter=0;
2741 //-- exclude the low intensity tracks
2743 for (Int_t index=0; index<indexMax; ++index){
2745 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2746 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2747 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2749 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2751 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2756 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2758 if (itrack>335) indexMax=557568/2;
2759 for (Int_t index=0; index<indexMax; ++index){
2760 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2761 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2762 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2763 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2765 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2766 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2767 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2768 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2770 //cut if no track info available
2771 if (iltr->GetBundle()==0) continue;
2772 if (iR<133||mR<133) continue;
2773 if(TMath::Abs(mltr->fVecP2->GetMatrixArray()[index])>minres[i]) continue;
2775 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2776 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2778 //Double_t xxx[4] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35};
2779 Double_t xxx[3] = {ldrift,iGy*ldrift/(zlength*250.), 250.-mR};
2781 if (iltr->GetSide()==0){
2782 fdriftA.AddPoint(xxx,mdrift,1);
2784 fdriftC.AddPoint(xxx,mdrift,1);
2786 // Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->fVecSec->GetMatrixArray()[index]>35, iltr->GetSide()};
2787 Double_t xxx2[4] = { ldrift,iGy*ldrift/(zlength*250.), 250.-mR, iltr->GetSide()};
2788 fdriftAC.AddPoint(xxx2,mdrift,1);
2791 }//end laser track loop
2801 fdriftA.GetParameters(fitA);
2802 fdriftC.GetParameters(fitC);
2803 fdriftAC.GetParameters(fitAC);
2805 //Parameters: 0 linear offset
2806 // 1 mean drift velocity correction factor
2807 // 2 relative global y gradient
2808 // 3 radial deformation
2809 // 4 IROC/OROC offset
2811 // FindResiduals(arrMeasured,arrIdeal,fitA,fitC);
2813 for (Int_t itrack=0; itrack<338; ++itrack){
2814 AliTPCLaserTrack *iltr=(AliTPCLaserTrack*)arrIdeal->UncheckedAt(itrack);
2815 AliTPCLaserTrack *mltr=(AliTPCLaserTrack*)arrMeasured->UncheckedAt(itrack);
2817 //-- Exclude the tracks which has the biggest inclanation angle
2818 if ((itrack%7==0||itrack%7==6)&&itrack<336) continue;
2819 Int_t clustercounter=0;
2822 //-- exclude the low intensity tracks
2824 for (Int_t index=0; index<indexMax; ++index){
2825 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2826 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2827 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2828 if (TMath::Abs(mGz)<1e-20 && TMath::Abs(mGy)<1e-20 && TMath::Abs(mGx)<1e-20) clustercounter++;
2830 if (clustercounter>130&&itrack<336) continue; // don't accept tracks with <= 159-130=29 clusters
2834 Double_t zlength = (iltr->GetSide()==0)? fParam->GetZLength(36): fParam->GetZLength(71);
2836 if (itrack>335) indexMax=557568/2;
2837 for (Int_t index=0; index<indexMax; ++index){
2838 Double_t iGx=iltr->fVecGX->GetMatrixArray()[index];
2839 Double_t iGy=iltr->fVecGY->GetMatrixArray()[index];
2840 Double_t iGz=iltr->fVecGZ->GetMatrixArray()[index];
2841 Double_t iR=TMath::Sqrt(iGx*iGx+iGy*iGy);
2843 Double_t mGx=mltr->fVecGX->GetMatrixArray()[index];
2844 Double_t mGy=mltr->fVecGY->GetMatrixArray()[index];
2845 Double_t mGz=mltr->fVecGZ->GetMatrixArray()[index];
2846 Double_t mR=TMath::Sqrt(mGx*mGx+mGy*mGy);
2848 //cut if no track info available
2849 if (iR<60||mR<60) continue;
2851 Double_t ldrift = (iltr->GetSide()==0)?zlength-iGz:iGz+zlength;
2852 Double_t mdrift = (iltr->GetSide()==0)?zlength-mGz:mGz+zlength;
2855 if (iltr->GetSide()==1) v=&fitC;
2856 // 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);
2857 Double_t iCorr=(*v)[0]+(*v)[1]*ldrift+(*v)[2]*iGy*ldrift/(zlength*250.)+(*v)[3]*(250.-mR);
2859 mltr->fVecP2->GetMatrixArray()[index]=mdrift-iCorr;
2868 //set statistics values
2870 npointsA= fdriftA.GetNpoints();
2871 if (npointsA>0) chi2A = fdriftA.GetChisquare()/fdriftA.GetNpoints();
2875 npointsC= fdriftC.GetNpoints();
2876 if (npointsC>0) chi2C = fdriftC.GetChisquare()/fdriftC.GetNpoints();
2880 npointsAC= fdriftAC.GetNpoints();
2881 if (npointsAC>0) chi2AC = fdriftAC.GetChisquare()/fdriftAC.GetNpoints();
2885 if (fStreamLevel>2){
2886 //laser track information
2887 (*GetDebugStreamer()) << "DriftV" <<
2889 "run=" << fRunNumber <<
2890 "time=" << timestamp <<
2891 "fitA.=" << &fitA <<
2892 "fitC.=" << &fitC <<
2893 "fitAC.=" << &fitAC <<
2903 //Parameters: 0 linear offset (global)
2904 // 1 mean drift velocity correction factor
2905 // 2 relative global y gradient
2906 // 3 radial deformation
2907 // 4 IROC/OROC offset
2908 // 5 linear offset relative A-C
2911 TGraphErrors *grA[7];
2912 TGraphErrors *grC[7];
2913 TGraphErrors *grAC[8];
2914 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_");
2915 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_");
2917 TObjArray *arrNames=names.Tokenize(";");
2918 TObjArray *arrNamesAC=namesAC.Tokenize(";");
2921 for (Int_t i=0; i<7; ++i){
2922 TString grName=arrNames->UncheckedAt(i)->GetName();
2924 grA[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2926 grA[i]=new TGraphErrors;
2927 grA[i]->SetName(grName.Data());
2928 grA[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2929 fArrFitGraphs->Add(grA[i]);
2931 // Int_t ipoint=grA[i]->GetN();
2933 grA[i]->SetPoint(ipoint,timestamp,fitA(i));
2934 grA[i]->SetPointError(ipoint,60,.0001);
2935 if (i<4) grA[i]->SetPointError(ipoint,60,fdriftA.GetCovarianceMatrixElement(i,i));
2939 for (Int_t i=0; i<7; ++i){
2940 TString grName=arrNames->UncheckedAt(i)->GetName();
2942 grC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2944 grC[i]=new TGraphErrors;
2945 grC[i]->SetName(grName.Data());
2946 grC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2947 fArrFitGraphs->Add(grC[i]);
2949 // Int_t ipoint=grC[i]->GetN();
2951 grC[i]->SetPoint(ipoint,timestamp,fitC(i));
2952 grC[i]->SetPointError(ipoint,60,.0001);
2953 if (i<4) grC[i]->SetPointError(ipoint,60,fdriftC.GetCovarianceMatrixElement(i,i));
2957 for (Int_t i=0; i<8; ++i){
2958 TString grName=arrNamesAC->UncheckedAt(i)->GetName();
2960 grAC[i]=(TGraphErrors*)fArrFitGraphs->FindObject(grName.Data());
2962 grAC[i]=new TGraphErrors;
2963 grAC[i]->SetName(grName.Data());
2964 grAC[i]->SetTitle(grName.ReplaceAll("_"," ").Data());
2965 fArrFitGraphs->Add(grAC[i]);
2967 // Int_t ipoint=grAC[i]->GetN();
2969 grAC[i]->SetPoint(ipoint,timestamp,fitAC(i));
2970 grAC[i]->SetPointError(ipoint,60,.0001);
2971 if (i<5) grAC[i]->SetPointError(ipoint,60,fdriftAC.GetCovarianceMatrixElement(i,i));
2974 if (fDebugLevel>10){
2975 printf("A side fit parameters:\n");
2977 printf("\nC side fit parameters:\n");
2979 printf("\nAC side fit parameters:\n");
2986 //_____________________________________________________________________
2987 Double_t AliTPCCalibCE::SetBurstHnDrift()
2990 // Create a new THnSparse for the current burst
2991 // return the time of the current burst
2994 for(i=0; i<fTimeBursts.GetNrows(); ++i){
2995 if(fTimeBursts.GetMatrixArray()[i]<1.e-20) break;
2996 if(TMath::Abs(fTimeBursts.GetMatrixArray()[i]-fTimeStamp)<300){
2997 fHnDrift=(THnSparseI*)fArrHnDrift.UncheckedAt(i);
2998 return fTimeBursts(i);
3003 fArrHnDrift.AddAt(fHnDrift,i);
3004 fTimeBursts.GetMatrixArray()[i]=fTimeStamp;
3013 //_____________________________________________________________________
3014 void AliTPCCalibCE::DumpToFile(const Char_t *filename, const Char_t *dir, Bool_t /*append*/)
3017 // Write class to file
3018 // option can be specified in the dir option:
3020 // name=<objname>: the name of the calibration object in file will be <objname>
3021 // type=<type>: the saving type:
3022 // 0 - write the complte object
3023 // 1 - Store the histogram arrays separately to make the streamed object smaller, Analyse to be called
3024 // 2 - like 2, but in addition delete objects that will most probably not be used for calibration
3025 // 3 - store only calibration output, don't store the reference histograms
3026 // and THnSparse (requires Analyse called before)
3028 // NOTE: to read the object back, the ReadFromFile function should be used
3032 TString objName=GetName();
3036 TObjArray *arr=sDir.Tokenize(",");
3039 while ( (s=(TObjString*)next()) ){
3040 TString optString=s->GetString();
3041 optString.Remove(TString::kBoth,' ');
3042 if (optString.BeginsWith("name=")){
3043 objName=optString.ReplaceAll("name=","");
3045 if (optString.BeginsWith("type=")){
3046 optString.ReplaceAll("type=","");
3047 type=optString.Atoi();
3052 // only for the new algorithm
3054 ce.fArrFitGraphs=fArrFitGraphs;
3055 ce.fNevents=fNevents;
3056 ce.fTimeBursts.ResizeTo(fTimeBursts.GetNrows());
3057 ce.fTimeBursts=fTimeBursts;
3058 ce.fProcessNew=kTRUE;
3059 TFile f(filename,"recreate");
3060 ce.Write(objName.Data());
3061 fArrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3063 ce.fArrFitGraphs=0x0;
3068 if (type==1||type==2) {
3069 //delete most probably not needed stuff
3070 //This requires Analyse to be called after reading the object from file
3071 fCalRocArrayT0.Delete();
3072 fCalRocArrayT0Err.Delete();
3073 fCalRocArrayQ.Delete();
3074 fCalRocArrayRMS.Delete();
3075 fCalRocArrayOutliers.Delete();
3077 if (type==2||type==3){
3078 fParamArrayEventPol1.Delete();
3079 fParamArrayEventPol2.Delete();
3082 TObjArray histoQArray(72);
3083 TObjArray histoT0Array(72);
3084 TObjArray histoRMSArray(72);
3085 TObjArray arrHnDrift(fArrHnDrift.GetEntries());
3087 //save all large 2D histograms in separte pointers
3088 //to have a smaller memory print when saving the object
3089 if (type==1||type==2||type==3){
3090 for (Int_t i=0; i<72; ++i){
3091 histoQArray.AddAt(fHistoQArray.UncheckedAt(i),i);
3092 histoT0Array.AddAt(fHistoT0Array.UncheckedAt(i),i);
3093 histoRMSArray.AddAt(fHistoRMSArray.UncheckedAt(i),i);
3095 fHistoQArray.SetOwner(kFALSE);
3096 fHistoT0Array.SetOwner(kFALSE);
3097 fHistoRMSArray.SetOwner(kFALSE);
3098 fHistoQArray.Clear();
3099 fHistoT0Array.Clear();
3100 fHistoRMSArray.Clear();
3102 for (Int_t i=0;i<fArrHnDrift.GetEntries();++i){
3103 arrHnDrift.AddAt(fArrHnDrift.UncheckedAt(i),i);
3105 fArrHnDrift.SetOwner(kFALSE);
3106 fArrHnDrift.Clear();
3110 TDirectory *backup = gDirectory;
3112 TFile f(filename,"recreate");
3113 Write(objName.Data());
3114 if (type==1||type==2) {
3115 histoQArray.Write("histoQArray",TObject::kSingleKey);
3116 histoT0Array.Write("histoT0Array",TObject::kSingleKey);
3117 histoRMSArray.Write("histoRMSArray",TObject::kSingleKey);
3118 arrHnDrift.Write("arrHnDrift",TObject::kSingleKey);
3124 //move histograms back to the object
3125 if (type==1||type==2){
3126 for (Int_t i=0; i<72; ++i){
3127 fHistoQArray.AddAt(histoQArray.UncheckedAt(i),i);
3128 fHistoT0Array.AddAt(histoT0Array.UncheckedAt(i),i);
3129 fHistoRMSArray.AddAt(histoRMSArray.UncheckedAt(i),i);
3131 fHistoQArray.SetOwner(kTRUE);
3132 fHistoT0Array.SetOwner(kTRUE);
3133 fHistoRMSArray.SetOwner(kTRUE);
3135 for (Int_t i=0;i<arrHnDrift.GetEntries();++i){
3136 fArrHnDrift.AddAt(arrHnDrift.UncheckedAt(i),i);
3138 fArrHnDrift.SetOwner(kTRUE);
3141 if ( backup ) backup->cd();
3143 //_____________________________________________________________________
3144 AliTPCCalibCE* AliTPCCalibCE::ReadFromFile(const Char_t *filename)
3147 // Read object from file
3148 // Handle properly if the histogram arrays were stored separately
3149 // call Analyse to make sure to have the calibration relevant information in the object
3153 if (!f.IsOpen() || f.IsZombie() ) return 0x0;
3154 TList *l=f.GetListOfKeys();
3159 AliTPCCalibCE *ce=0x0;
3160 TObjArray *histoQArray=0x0;
3161 TObjArray *histoT0Array=0x0;
3162 TObjArray *histoRMSArray=0x0;
3163 TObjArray *arrHnDrift=0x0;
3165 while ( (key=(TKey*)next()) ){
3167 if ( o->IsA()==AliTPCCalibCE::Class() ){
3168 ce=(AliTPCCalibCE*)o;
3169 } else if ( o->IsA()==TObjArray::Class() ){
3170 TString name=key->GetName();
3171 if ( name=="histoQArray") histoQArray=(TObjArray*)o;
3172 if ( name=="histoT0Array") histoT0Array=(TObjArray*)o;
3173 if ( name=="histoRMSArray") histoRMSArray=(TObjArray*)o;
3174 if ( name=="arrHnDrift") arrHnDrift=(TObjArray*)o;
3179 //move histograms back to the object
3182 for (Int_t i=0; i<72; ++i){
3183 hist=(TH1*)histoQArray->UncheckedAt(i);
3184 if (hist) hist->SetDirectory(0x0);
3185 ce->fHistoQArray.AddAt(hist,i);
3187 ce->fHistoQArray.SetOwner(kTRUE);
3191 for (Int_t i=0; i<72; ++i){
3192 hist=(TH1*)histoT0Array->UncheckedAt(i);
3193 if (hist) hist->SetDirectory(0x0);
3194 ce->fHistoT0Array.AddAt(hist,i);
3196 ce->fHistoT0Array.SetOwner(kTRUE);
3200 for (Int_t i=0; i<72; ++i){
3201 hist=(TH1*)histoRMSArray->UncheckedAt(i);
3202 if (hist) hist->SetDirectory(0x0);
3203 ce->fHistoRMSArray.AddAt(hist,i);
3205 ce->fHistoRMSArray.SetOwner(kTRUE);
3209 for (Int_t i=0; i<arrHnDrift->GetEntries(); ++i){
3210 THnSparseI *hSparse=(THnSparseI*)arrHnDrift->UncheckedAt(i);
3211 ce->fArrHnDrift.AddAt(hSparse,i);