]> git.uio.no Git - u/mrichter/AliRoot.git/blame - T0/AliT0Reconstructor.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
new pedestal treatment
[u/mrichter/AliRoot.git] / T0 / AliT0Reconstructor.cxx
CommitLineData
dc7ca31d 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16/* $Id$ */
72e48d95 17/*********************************************************************
18 * T0 reconstruction and filling ESD
19 * - reconstruct mean time (interation time)
20 * - vertex position
21 * - multiplicity
22 ********************************************************************/
dc7ca31d 23
af885e0f 24#include <AliESDEvent.h>
dc7ca31d 25#include "AliLog.h"
dc7ca31d 26#include "AliT0RecPoint.h"
27#include "AliRawReader.h"
28#include "AliT0RawReader.h"
dc7ca31d 29#include "AliT0digit.h"
30#include "AliT0Reconstructor.h"
31#include "AliT0Parameters.h"
c41ceaac 32#include "AliT0Calibrator.h"
58bd3a16 33#include "AliESDfriend.h"
b0e13b29 34#include "AliESDTZERO.h"
73df58ab 35#include "AliESDTZEROfriend.h"
8f8d0732 36#include "AliLog.h"
85f61e3b 37#include "AliCDBEntry.h"
38#include "AliCDBManager.h"
39#include "AliCTPTimeParams.h"
40#include "AliLHCClockPhase.h"
669dc07f 41#include "AliT0CalibSeasonTimeShift.h"
4cbe597e 42#include "AliESDRun.h"
dc7ca31d 43
44#include <TArrayI.h>
45#include <TGraph.h>
aad72f45 46#include <TMath.h>
b09247a2 47#include <Riostream.h>
dc7ca31d 48
49ClassImp(AliT0Reconstructor)
50
c41ceaac 51 AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
f16935f7 52 fdZonA(0),
53 fdZonC(0),
54 fZposition(0),
55 fParam(NULL),
2e6a5ee0 56 fAmpLEDrec(),
c883fdf2 57 fQTC(0),
58 fAmpLED(0),
58bd3a16 59 fCalib(),
60 fLatencyHPTDC(9000),
61 fLatencyL1(0),
62 fLatencyL1A(0),
73df58ab 63 fLatencyL1C(0),
85f61e3b 64 fGRPdelays(0),
669dc07f 65 fTimeMeanShift(0x0),
66 fTimeSigmaShift(0x0),
b0e13b29 67 fESDTZEROfriend(NULL),
68 fESDTZERO(NULL)
58bd3a16 69
e0bba6cc 70{
612737bb 71 for (Int_t i=0; i<24; i++) fTime0vertex[i] =0;
4cbe597e 72
72e48d95 73 //constructor
85f61e3b 74 AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/CTP/CTPtiming");
75 if (!entry) AliFatal("CTP timing parameters are not found in OCDB !");
76 AliCTPTimeParams *ctpParams = (AliCTPTimeParams*)entry->GetObject();
77 Float_t l1Delay = (Float_t)ctpParams->GetDelayL1L0()*25.0;
78
79 AliCDBEntry *entry1 = AliCDBManager::Instance()->Get("GRP/CTP/TimeAlign");
80 if (!entry1) AliFatal("CTP time-alignment is not found in OCDB !");
81 AliCTPTimeParams *ctpTimeAlign = (AliCTPTimeParams*)entry1->GetObject();
82 l1Delay += ((Float_t)ctpTimeAlign->GetDelayL1L0()*25.0);
83
84 AliCDBEntry *entry4 = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase");
85 if (!entry4) AliFatal("LHC clock-phase shift is not found in OCDB !");
86 AliLHCClockPhase *phase = (AliLHCClockPhase*)entry4->GetObject();
87
85f61e3b 88 fGRPdelays = l1Delay - phase->GetMeanPhase();
72e48d95 89
669dc07f 90 AliCDBEntry *entry5 = AliCDBManager::Instance()->Get("T0/Calib/TimeAdjust");
91 if (entry5) {
92 AliT0CalibSeasonTimeShift *timeshift = (AliT0CalibSeasonTimeShift*)entry5->GetObject();
93 fTimeMeanShift = timeshift->GetT0Means();
94 fTimeSigmaShift = timeshift->GetT0Sigmas();
4cbe597e 95 }
669dc07f 96 else
97 AliWarning("Time Adjust is not found in OCDB !");
612737bb 98
74adb36a 99 fParam = AliT0Parameters::Instance();
100 fParam->Init();
c883fdf2 101
74adb36a 102 for (Int_t i=0; i<24; i++){
2e6a5ee0 103 TGraph* gr = fParam ->GetAmpLEDRec(i);
29ed1d0e 104 if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ;
c883fdf2 105 TGraph* gr1 = fParam ->GetAmpLED(i);
106 if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ;
107 TGraph* gr2 = fParam ->GetQTC(i);
539b9cb9 108 if (gr2) fQTC.AddAtAndExpand(gr2,i) ;
612737bb 109 fTime0vertex[i] = fParam->GetCFD(i);
b0e13b29 110 AliDebug(2,Form("OCDB mean CFD time %i %f \n",i, fTime0vertex[i]));
612737bb 111 }
58bd3a16 112 fLatencyL1 = fParam->GetLatencyL1();
612737bb 113 fLatencyL1A = fParam->GetLatencyL1A();
58bd3a16 114 fLatencyL1C = fParam->GetLatencyL1C();
115 fLatencyHPTDC = fParam->GetLatencyHPTDC();
9d026202 116 AliDebug(2,Form(" LatencyL1 %f latencyL1A %f latencyL1C %f latencyHPTDC %f \n",fLatencyL1, fLatencyL1A, fLatencyL1C, fLatencyHPTDC));
612737bb 117
118 for (Int_t i=0; i<24; i++) {
59fe6376 119 if( fTime0vertex[i] < 500 || fTime0vertex[i] > 50000) fTime0vertex[i] =( 1000.*fLatencyHPTDC - 1000.*fLatencyL1 + 1000.*fGRPdelays)/24.4;
92f6bfd3 120 // printf(" calulated mean %i %f \n",i, fTime0vertex[i]);
612737bb 121 }
8f620945 122 //here real Z position
123 fdZonC = TMath::Abs(fParam->GetZPosition("T0/C/PMT1"));
124 fdZonA = TMath::Abs(fParam->GetZPosition("T0/A/PMT15"));
b7fd8e80 125 // printf(" !!!! fdZonC %f fdZonA %f \n", fdZonC/29.9792458, fdZonA/29.9792458);
12e9daf9 126 fCalib = new AliT0Calibrator();
73df58ab 127 fESDTZEROfriend = new AliESDTZEROfriend();
b0e13b29 128 fESDTZERO = new AliESDTZERO();
12e9daf9 129
b0e13b29 130
dc7ca31d 131}
c41ceaac 132
133//_____________________________________________________________________________
dc7ca31d 134void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
135{
94c27e4f 136 // T0 digits reconstruction
ce50812a 137 Int_t refAmp = 0 ; /*Int_t (GetRecoParam()->GetRefAmp());*/
776de217 138
c41ceaac 139 TArrayI * timeCFD = new TArrayI(24);
140 TArrayI * timeLED = new TArrayI(24);
141 TArrayI * chargeQT0 = new TArrayI(24);
142 TArrayI * chargeQT1 = new TArrayI(24);
74adb36a 143
d0bcd1fb 144
ce50812a 145 Float_t c = 29.9792458; // cm/ns
8955c6b4 146 Float_t channelWidth = fParam->GetChannelWidth() ;
eadc5bd4 147 Double32_t vertex = 9999999, meanVertex = 0 ;
adf36b9d 148 Double32_t timeDiff=999999, meanTime=999999, timeclock=999999;
ce50812a 149
94c27e4f 150
dc7ca31d 151 AliDebug(1,Form("Start DIGITS reconstruction "));
94c27e4f 152
ce50812a 153 Float_t lowAmpThreshold = GetRecoParam()->GetAmpLowThreshold();
154 Float_t highAmpThreshold = GetRecoParam()->GetAmpHighThreshold();
155 printf( "AliT0Reconstructor::Reconstruct::: RecoParam amplitude %f %f \n",lowAmpThreshold, highAmpThreshold);
156
157 Double32_t besttimeA=9999999; Double32_t besttimeA_best=9999999;
158 Double32_t besttimeC=9999999; Double32_t besttimeC_best=9999999;
159 Int_t timeDelayCFD[24];
160 Int_t badpmt[24];
161 //Bad channel
162 for (Int_t i=0; i<24; i++) {
163 badpmt[i] = GetRecoParam() -> GetBadChannels(i);
164 timeDelayCFD[i] = Int_t (fParam->GetTimeDelayCFD(i));
165 }
166 fCalib->SetEq(0);
d0bcd1fb 167 TBranch *brDigits=digitsTree->GetBranch("T0");
e0bba6cc 168 AliT0digit *fDigits = new AliT0digit() ;
dc7ca31d 169 if (brDigits) {
170 brDigits->SetAddress(&fDigits);
171 }else{
f16935f7 172 AliError(Form("EXEC Branch T0 digits not found"));
ce50812a 173 return;
dc7ca31d 174 }
e0bba6cc 175
c41ceaac 176 digitsTree->GetEvent(0);
177 digitsTree->GetEntry(0);
178 brDigits->GetEntry(0);
179 fDigits->GetTimeCFD(*timeCFD);
180 fDigits->GetTimeLED(*timeLED);
181 fDigits->GetQT0(*chargeQT0);
182 fDigits->GetQT1(*chargeQT1);
446d6ec4 183 Int_t onlineMean = fDigits->MeanTime();
ce50812a 184
adf36b9d 185 Bool_t tr[5];
186 for (Int_t i=0; i<5; i++) tr[i]=false;
c41ceaac 187
dc7ca31d 188
7e08771f 189 AliT0RecPoint frecpoints;
190 AliT0RecPoint * pfrecpoints = &frecpoints;
191 clustersTree->Branch( "T0", "AliT0RecPoint" ,&pfrecpoints);
94c27e4f 192
195e1353 193 Float_t time[24], adc[24], adcmip[24];
dc7ca31d 194 for (Int_t ipmt=0; ipmt<24; ipmt++) {
ce50812a 195 if(timeCFD->At(ipmt)>0 ) {
b7fd8e80 196 Float_t timefull = 0.001*( timeCFD->At(ipmt) - 513 - timeDelayCFD[ipmt]) * channelWidth;
eadc5bd4 197 frecpoints.SetTimeFull(ipmt, 0 ,timefull) ;
ce50812a 198 if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)
d0bcd1fb 199 adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
a7027400 200 else
d0bcd1fb 201 adc[ipmt] = 0;
202
ce50812a 203 time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]), timeCFD->At(ipmt)) ;
b7fd8e80 204 time[ipmt] = time[ipmt] - 513;
d0bcd1fb 205 Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
8f620945 206 // time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD->At(ipmt) ) ;
eadc5bd4 207 AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
ce50812a 208 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
209
644a358e 210 Double_t ampMip = 0;
211 TGraph* ampGraph = (TGraph*)fAmpLED.At(ipmt);
212 if (ampGraph) ampMip = ampGraph->Eval(sl);
213 Double_t qtMip = 0;
214 TGraph* qtGraph = (TGraph*)fQTC.At(ipmt);
215 if (qtGraph) qtMip = qtGraph->Eval(adc[ipmt]);
669dc07f 216 AliDebug(5,Form(" Amlitude in MIPS LED %f , QTC %f in channels %f\n ",ampMip,qtMip, adc[ipmt]));
7e08771f 217 frecpoints.SetTime(ipmt, Float_t(time[ipmt]) );
218 frecpoints.SetAmpLED(ipmt, Float_t( ampMip));
219 frecpoints.SetAmp(ipmt, Float_t(qtMip));
195e1353 220 adcmip[ipmt]=qtMip;
d0bcd1fb 221
dc7ca31d 222 }
223 else {
eadc5bd4 224 time[ipmt] = -99999;
dc7ca31d 225 adc[ipmt] = 0;
b0ab3f59 226 adcmip[ipmt] = 0;
ce50812a 227
dc7ca31d 228 }
229 }
94c27e4f 230
dc7ca31d 231 for (Int_t ipmt=0; ipmt<12; ipmt++){
ce50812a 232 if(time[ipmt] !=0 && time[ipmt] > -9000
233 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold )
234 {
235 if(time[ipmt]<besttimeC) besttimeC=time[ipmt]; //timeC
236 if(TMath::Abs(time[ipmt])<TMath::Abs(besttimeC_best))
237 besttimeC_best=time[ipmt]; //timeC
dc7ca31d 238 }
dc7ca31d 239 }
ce50812a 240 for ( Int_t ipmt=12; ipmt<24; ipmt++)
241 {
242 if(time[ipmt] != 0 && time[ipmt] > -9000
243 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold)
244 {
245 if(time[ipmt]<besttimeA) besttimeA=time[ipmt];
246 if(TMath::Abs(time[ipmt] ) < TMath::Abs(besttimeA_best))
247 besttimeA_best=time[ipmt]; //timeA
248 }
dc7ca31d 249 }
ce50812a 250
251 if( besttimeA < 999999 && besttimeA!=0) {
eadc5bd4 252 frecpoints.SetTimeBestA((besttimeA_best * channelWidth - fdZonA/c) );
b7fd8e80 253 frecpoints.SetTime1stA((besttimeA * channelWidth - fdZonA/c ) );
adf36b9d 254 tr[1]=true;
255 }
ce50812a 256
257 if( besttimeC < 999999 && besttimeC!=0) {
eadc5bd4 258 frecpoints.SetTimeBestC((besttimeC_best * channelWidth - fdZonC/c) );
b7fd8e80 259 frecpoints.SetTime1stC((besttimeC * channelWidth - fdZonC/c ) );
adf36b9d 260 tr[2]=true;
261 }
ce50812a 262
263 AliDebug(5,Form("1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f ",
264 besttimeA, besttimeA_best,
265 besttimeC, besttimeC_best) );
266
adf36b9d 267 if(besttimeA <999999 && besttimeC < 999999 ){
9b83615d 268 // timeDiff = (besttimeC - besttimeA)*channelWidth;
269 timeDiff = (besttimeA - besttimeC)*channelWidth;
ce50812a 270 meanTime = channelWidth * (besttimeA_best + besttimeC_best)/2. ;
271 timeclock = channelWidth * (besttimeA + besttimeC)/2. ;
eadc5bd4 272 vertex = meanVertex - 0.001* c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
adf36b9d 273 tr[0]=true;
274 }
7e08771f 275 frecpoints.SetVertex(vertex);
eadc5bd4 276 frecpoints.SetMeanTime(meanTime );
277 frecpoints.SetT0clock(timeclock );
7e08771f 278 frecpoints.SetT0Trig(tr);
ce50812a 279
eadc5bd4 280 AliDebug(5,Form("fRecPoints::: 1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f vertex %f",
ce50812a 281 frecpoints.Get1stTimeA(), frecpoints.GetBestTimeA(),
282 frecpoints.Get1stTimeC(), frecpoints.GetBestTimeC(),
eadc5bd4 283 vertex ) );
ce50812a 284
669dc07f 285 AliDebug(5,Form("T0 triggers %d %d %d %d %d",tr[0],tr[1],tr[2],tr[3],tr[4]));
ce50812a 286
adf36b9d 287 //online mean
7e08771f 288 frecpoints.SetOnlineMean(Int_t(onlineMean));
b5a9f753 289 AliDebug(10,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm online mean %i timeclock %f ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
adf36b9d 290
830172c8 291 clustersTree->Fill();
ce50812a 292
bd375212 293 delete timeCFD;
294 delete timeLED;
295 delete chargeQT0;
296 delete chargeQT1;
dc7ca31d 297}
298
299
c41ceaac 300//_______________________________________________________________________
301
302void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
303{
94c27e4f 304 // T0 raw ->
539b9cb9 305 //
c8d939c8 306
b0e13b29 307 Float_t meanOrA = fTime0vertex[0] + 587;
308 Float_t meanOrC = fTime0vertex[0] + 678;
309 Float_t meanTVDC = fTime0vertex[0] + 2564;
82540e64 310 Int_t timeDelayCFD[24];
311
d3e04608 312
612737bb 313 Int_t badpmt[24];
38cbfa7c 314 //Bad channel
58432641 315 for (Int_t i=0; i<24; i++) {
316 badpmt[i] = GetRecoParam() -> GetBadChannels(i);
82540e64 317 timeDelayCFD[i] = Int_t (fParam->GetTimeDelayCFD(i));
58432641 318 }
ce50812a 319 Int_t equalize = GetRecoParam() -> GetEq();
c8d939c8 320 // printf( "AliT0Reconstructor::Reconstruct::: RecoParam %i \n",equalize);
ce50812a 321 fCalib->SetEq(equalize);
d3e04608 322 Int_t low[500], high[500];
94b4f7e2 323 Float_t timefull=-9999;;
324 Float_t tvdc = -9999; Float_t ora = -9999; Float_t orc = -9999;
58bd3a16 325
e8ed1cd0 326 Int_t allData[110][5];
2e6a5ee0 327
e8ed1cd0 328 Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
ce50812a 329 Float_t time2zero[24];
612737bb 330 Double32_t timeDiff, meanTime, timeclock;
331 timeDiff = meanTime = timeclock = 9999999;
adf36b9d 332 Float_t c = 29.9792458; // cm/ns
333 Double32_t vertex = 9999999;
776de217 334 Int_t onlineMean=0;
9480f05f 335 Float_t meanVertex = 0;
b7fd8e80 336 Int_t pedestal[24];
612737bb 337 for (Int_t i0=0; i0<24; i0++) {
b7fd8e80 338 low[i0] = Int_t(fTime0vertex[i0]) - 50;
339 high[i0] = Int_t(fTime0vertex[i0]) + 50;
340 time2zero[i0] = 999999;
341 pedestal[i0]=GetRecoParam()->GetLow(i0+100);
342 // printf("pmt %i pedestal %i\n", i0,pedestal[i0]);
612737bb 343 }
344
669dc07f 345 for (Int_t i0=0; i0<110; i0++)
b0e13b29 346 for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
612737bb 347
348 Float_t lowAmpThreshold = GetRecoParam()->GetAmpLowThreshold();
349 Float_t highAmpThreshold = GetRecoParam()->GetAmpHighThreshold();
b7fd8e80 350 // printf( "AliT0Reconstructor::Reconstruct::: RecoParam amplitude %f %f \n",lowAmpThreshold, highAmpThreshold);
ce50812a 351
352 Double32_t besttimeA=9999999; Double32_t besttimeA_best=9999999;
353 Double32_t besttimeC=9999999; Double32_t besttimeC_best=9999999;
354
355 Float_t channelWidth = fParam->GetChannelWidth() ;
b0e13b29 356
7e08771f 357 AliT0RecPoint frecpoints;
358 AliT0RecPoint * pfrecpoints = &frecpoints;
bce12dc5 359
7e08771f 360 recTree->Branch( "T0", "AliT0RecPoint" ,&pfrecpoints);
2e6a5ee0 361
bce12dc5 362 AliDebug(10," before read data ");
363 AliT0RawReader myrawreader(rawReader);
776de217 364
365 UInt_t type =rawReader->GetType();
366
bce12dc5 367 if (!myrawreader.Next())
368 AliDebug(1,Form(" no raw data found!!"));
369 else
370 {
38cbfa7c 371 for (Int_t i=0; i<24; i++)
372 {
373 timeCFD[i]=0; timeLED[i]=0; chargeQT0[i]=0; chargeQT1[i]=0;
374 }
375 Int_t fBCID=Int_t (rawReader->GetBCID());
8f620945 376 Int_t trmbunch= myrawreader.GetTRMBunchID();
86fd0587 377 AliDebug(10,Form(" CDH BC ID %i, TRM BC ID %i \n", fBCID, trmbunch ));
378
38cbfa7c 379 if(type == 7 ) { //only physics
669dc07f 380 for (Int_t i=0; i<107; i++) {
bce12dc5 381 for (Int_t iHit=0; iHit<5; iHit++)
382 {
383 allData[i][iHit] = myrawreader.GetData(i,iHit);
384 }
8f620945 385 }
8f620945 386
85f61e3b 387 for (Int_t in=0; in<12; in++)
388 {
389 for (Int_t iHit=0; iHit<5; iHit++)
38cbfa7c 390 {
612737bb 391 if(allData[in+1][iHit] > low[in] &&
392 allData[in+1][iHit] < high[in])
38cbfa7c 393 {
85f61e3b 394 timeCFD[in] = allData[in+1][iHit] ;
395 break;
38cbfa7c 396 }
85f61e3b 397 }
398 for (Int_t iHit=0; iHit<5; iHit++)
399 {
612737bb 400 if(allData[in+1+56][iHit] > low[in] &&
401 allData[in+1+56][iHit] < high[in])
38cbfa7c 402 {
85f61e3b 403 timeCFD[in+12] = allData[in+56+1][iHit] ;
404 break;
38cbfa7c 405 }
38cbfa7c 406 }
612737bb 407 timeLED[in+12] = allData[in+68+1][0] ;
408 timeLED[in] = allData[in+12+1][0] ;
ce50812a 409 AliDebug(50, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
85f61e3b 410 in, timeCFD[in],timeCFD[in+12],timeLED[in],
612737bb 411 timeLED[in+12]));
38cbfa7c 412
8f620945 413 }
414
8f620945 415
85f61e3b 416 for (Int_t in=0; in<12; in++)
417 {
418 chargeQT0[in]=allData[2*in+25][0];
419 chargeQT1[in]=allData[2*in+26][0];
612737bb 420 AliDebug(25, Form(" readed Raw %i %i %i",
85f61e3b 421 in, chargeQT0[in],chargeQT1[in]));
422 }
423 for (Int_t in=12; in<24; in++)
424 {
425 chargeQT0[in]=allData[2*in+57][0];
426 chargeQT1[in]=allData[2*in+58][0];
612737bb 427 AliDebug(25, Form(" readed Raw %i %i %i",
85f61e3b 428 in, chargeQT0[in],chargeQT1[in]));
85f61e3b 429 }
430
8f226345 431 onlineMean = allData[49][0];
432
195e1353 433 Double32_t time[24], adc[24], adcmip[24], noncalibtime[24];
8f620945 434 for (Int_t ipmt=0; ipmt<24; ipmt++) {
b7fd8e80 435 if(timeCFD[ipmt] > 0 &&( chargeQT0[ipmt] - chargeQT1[ipmt])>pedestal[ipmt] /* && badpmt[ipmt]==0*/ ) {
d0bcd1fb 436 //for physics data
b7fd8e80 437 if(( chargeQT0[ipmt] - chargeQT1[ipmt])>pedestal[ipmt]) {
541b42c4 438 adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
73df58ab 439 }
d0bcd1fb 440 else
441 adc[ipmt] = 0;
612737bb 442 // time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
c8d939c8 443 Int_t refAmp = Int_t (fTime0vertex[ipmt]);
444 time[ipmt] = fCalib-> WalkCorrection( refAmp, ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
d0bcd1fb 445 Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
8f620945 446 // time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
73df58ab 447 AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
d0bcd1fb 448 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
644a358e 449 Double_t ampMip = 0;
450 TGraph * ampGraph = (TGraph*)fAmpLED.At(ipmt);
451 if (ampGraph) ampMip = ampGraph->Eval(sl);
452 Double_t qtMip = 0;
453 TGraph * qtGraph = (TGraph*)fQTC.At(ipmt);
454 if (qtGraph) qtMip = qtGraph->Eval(adc[ipmt]);
345f03db 455 AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %f\n ",ampMip,qtMip, adc[ipmt]));
c8d939c8 456 if( qtMip>lowAmpThreshold && qtMip<highAmpThreshold )
457 {
458 if( equalize ==0 )
459 frecpoints.SetTime(ipmt, Float_t(time[ipmt]) );
ce50812a 460 else
461 frecpoints.SetTime(ipmt, Float_t(time[ipmt] + fTime0vertex[ipmt]) );
c8d939c8 462 // frecpoints.SetTime(ipmt, Float_t(time[ipmt] ) );
b7fd8e80 463
ce50812a 464 frecpoints.SetAmp(ipmt, Double32_t( qtMip));
465 adcmip[ipmt]=qtMip;
b7fd8e80 466 frecpoints.SetAmpLED(ipmt, Double32_t(ampMip));
467 }
ce50812a 468 noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
2e6a5ee0 469 }
470 else {
ce50812a 471 time[ipmt] = -9999;
2e6a5ee0 472 adc[ipmt] = 0;
b0ab3f59 473 adcmip[ipmt] = 0;
ce50812a 474 noncalibtime[ipmt] = -9999;
2e6a5ee0 475 }
476 }
ce50812a 477 fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;
478
2e6a5ee0 479 for (Int_t ipmt=0; ipmt<12; ipmt++){
ce50812a 480 if(time[ipmt] !=0 && time[ipmt] > -9000
481 /*&& badpmt[ipmt]==0 */
482 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold )
73df58ab 483 {
ce50812a 484 if(time[ipmt]<besttimeC) besttimeC=time[ipmt]; //timeC
485 if(TMath::Abs(time[ipmt])<TMath::Abs(besttimeC_best))
486 besttimeC_best=time[ipmt]; //timeC
2e6a5ee0 487 }
2e6a5ee0 488 }
73df58ab 489 for ( Int_t ipmt=12; ipmt<24; ipmt++)
490 {
ce50812a 491 if(time[ipmt] != 0 && time[ipmt] > -9000
492 /* && badpmt[ipmt]==0*/
b7fd8e80 493 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<40)
494 // && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold)
73df58ab 495 {
ce50812a 496 if(time[ipmt]<besttimeA) besttimeA=time[ipmt];
497 if(TMath::Abs(time[ipmt] ) < TMath::Abs(besttimeA_best))
498 besttimeA_best=time[ipmt]; //timeA
73df58ab 499 }
2e6a5ee0 500 }
ce50812a 501
502 if(besttimeA < 999999 && besttimeA!=0 ) {
503 if( equalize ==0 )
504 frecpoints.SetTime1stA((besttimeA * channelWidth)- 1000.*fLatencyHPTDC + 1000.*fLatencyL1A - 1000.*fGRPdelays - fTimeMeanShift[1] );
505 else
506 {
507 frecpoints.SetTimeBestA((besttimeA_best * channelWidth ));
508 frecpoints.SetTime1stA((besttimeA * channelWidth - fTimeMeanShift[1]));
509 }
510 }
511 if( besttimeC < 999999 && besttimeC!=0) {
512 if( equalize ==0 )
513 frecpoints.SetTime1stC((besttimeC * channelWidth)- 1000.*fLatencyHPTDC +1000.*fLatencyL1C - 1000.*fGRPdelays - fTimeMeanShift[2]);
514 else
515 {
516 frecpoints.SetTimeBestC((besttimeC_best * channelWidth ));
517 frecpoints.SetTime1stC((besttimeC * channelWidth - fTimeMeanShift[2]));
518 }
519 }
520 AliDebug(5,Form("1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f ",
521 besttimeA, besttimeA_best,
522 besttimeC, besttimeC_best) );
523 AliDebug(5,Form("fRecPoints::: 1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f shiftA %f shiftC %f ",
524 frecpoints.Get1stTimeA(), frecpoints.GetBestTimeA(),
525 frecpoints.Get1stTimeC(), frecpoints.GetBestTimeC(),
526 fTimeMeanShift[1],fTimeMeanShift[2] ) );
527 if( besttimeC < 999999 && besttimeA < 999999) {
c8d939c8 528 if(equalize ==0 )
529 timeclock = ((besttimeC * channelWidth)- 1000.*fLatencyHPTDC +1000.*fLatencyL1 - 1000.*fGRPdelays - fTimeMeanShift[0]);
530 else
531 {
532 timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - fTimeMeanShift[0];
533 meanTime = channelWidth * Float_t(besttimeA_best + besttimeC_best )/2.;
534 }
612737bb 535 timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth ;
9480f05f 536 vertex = meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2;
ce50812a 537 }
538
776de217 539 } //if phys event
ce50812a 540 AliDebug(1,Form(" timeDiff %f #channel, meanTime %f #ps, TOFmean%f vertex %f cm meanVertex %f \n",timeDiff, meanTime,timeclock, vertex,meanVertex));
541 frecpoints.SetT0clock(timeclock);
542 frecpoints.SetVertex(vertex);
543 frecpoints.SetMeanTime(meanTime);
544 frecpoints.SetOnlineMean(Int_t(onlineMean));
adf36b9d 545
ce50812a 546 // Set triggers
adf36b9d 547 Bool_t tr[5];
92f6bfd3 548 Int_t trchan[5] = {50,51,52,55,56};
549 Float_t lowtr[5] = {meanTVDC-700, meanOrA-700, meanOrC-700, meanOrC-1000, meanOrC-1000 };
550 Float_t hightr[5] = {meanTVDC+700, meanOrA+700, meanOrC+700, meanOrC+1000, meanOrC+1000};
ce50812a 551
82540e64 552 for (Int_t i=0; i<5; i++) tr[i] = false;
adf36b9d 553 for (Int_t itr=0; itr<5; itr++) {
b0e13b29 554 for (Int_t iHit=0; iHit<1; iHit++)
38cbfa7c 555 {
556 Int_t trr=trchan[itr];
92f6bfd3 557 if( allData[trr][iHit] > lowtr[itr] && allData[trr][iHit] < hightr[itr]) tr[itr]=true;
ce50812a 558
559 AliDebug(15,Form("Reconstruct ::: T0 triggers iHit %i tvdc %d orA %d orC %d centr %d semicentral %d",iHit, tr[0],tr[1],tr[2],tr[3],tr[4]));
b0e13b29 560 }
adf36b9d 561 }
7e08771f 562 frecpoints.SetT0Trig(tr);
b0e13b29 563
ce50812a 564 // all times with amplitude correction
565 Float_t timecent;
b0e13b29 566 for (Int_t iHit=0; iHit<5; iHit++)
567 {
ce50812a 568 timefull = timecent = -9999;
82540e64 569 tvdc = ora = orc = -9999;
b0e13b29 570 if(allData[50][iHit]>0)
571 tvdc = (Float_t(allData[50][iHit]) - meanTVDC) * channelWidth* 0.001;
572 if(allData[51][iHit]>0)
573 ora = (Float_t(allData[51][iHit]) - meanOrA) * channelWidth* 0.001;
574
575 if(allData[52][iHit]>0)
576 orc = (Float_t(allData[52][iHit]) - meanOrC) * channelWidth* 0.001;
577
7e08771f 578 frecpoints.SetOrC( iHit, orc);
579 frecpoints.SetOrA( iHit, ora);
580 frecpoints.SetTVDC( iHit, tvdc);
b0e13b29 581 for (Int_t i0=0; i0<12; i0++) {
ce50812a 582 if (equalize ==0 )
583 timecent = fTime0vertex[i0] + timeDelayCFD[i0];
584 else
585 timecent = fTime0vertex[i0];
b0e13b29 586 timefull = -9999;
587 if(allData[i0+1][iHit]>1)
ce50812a 588 timefull = (Float_t(allData[i0+1][iHit]) - timecent)* channelWidth* 0.001;
7e08771f 589 frecpoints.SetTimeFull(i0, iHit,timefull) ;
82540e64 590 // if(allData[i0+1][iHit]>1) printf("i0 %d iHit %d data %d fTime0vertex %f timefull %f \n",i0, iHit, allData[i0+1][iHit], fTime0vertex[i0], timefull);
b0e13b29 591
592 }
593
594 for (Int_t i0=12; i0<24; i0++) {
595 timefull = -9999;
ce50812a 596 if (equalize ==0 )
597 timecent = fTime0vertex[i0] + timeDelayCFD[i0];
598 else
599 timecent = fTime0vertex[i0];
94b4f7e2 600 if(allData[i0+45][iHit]>1) {
ce50812a 601 timefull = (Float_t(allData[i0+45][iHit]) - timecent)* channelWidth* 0.001;
94b4f7e2 602 }
82540e64 603 // if(allData[i0+45][iHit]>1) printf("i0 %d iHit %d data %d fTime0vertex %f timefull %f \n",i0, iHit, allData[i0+45][iHit], fTime0vertex[i0], timefull);
7e08771f 604 frecpoints.SetTimeFull(i0, iHit, timefull) ;
b0e13b29 605 }
606 }
607
608
669dc07f 609 //Set MPD
610 if(allData[53][0]>0 && allData[54][0])
7e08771f 611 frecpoints.SetMultA(allData[53][0]-allData[54][0]);
b0e13b29 612 if(allData[105][0]>0 && allData[106][0])
7e08771f 613 frecpoints.SetMultC(allData[105][0]-allData[106][0]);
b0e13b29 614
615
669dc07f 616 } // if (else )raw data
58bd3a16 617 recTree->Fill();
58bd3a16 618}
adf36b9d 619
620
ce50812a 621//____________________________________________________________
adf36b9d 622
623 void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
624 {
dc7ca31d 625
626 /***************************************************
627 Resonstruct digits to vertex position
628 ****************************************************/
629
dc7ca31d 630 AliDebug(1,Form("Start FillESD T0"));
b0ab3f59 631 if(!pESD) {
632 AliError("No ESD Event");
633 return;
634 }
4cbe597e 635 pESD ->SetT0spread(fTimeSigmaShift);
636
36cde487 637
58bd3a16 638 Float_t channelWidth = fParam->GetChannelWidth() ;
f7c2c2fc 639 Float_t c = 0.0299792458; // cm/ps
adf36b9d 640 Float_t currentVertex=0, shift=0;
291f31a1 641 Int_t ncont=-1;
adf36b9d 642 const AliESDVertex* vertex = pESD->GetPrimaryVertex();
643 if (!vertex) vertex = pESD->GetPrimaryVertexSPD();
644 if (!vertex) vertex = pESD->GetPrimaryVertexTPC();
645 if (!vertex) vertex = pESD->GetVertex();
646
647 if (vertex) {
648 AliDebug(2, Form("Got %s (%s) from ESD: %f",
649 vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
650 currentVertex = vertex->GetZ();
651
652 ncont = vertex->GetNContributors();
291f31a1 653 if(ncont>0 ) {
adf36b9d 654 shift = currentVertex/c;
adf36b9d 655 }
656 }
d76c31f4 657 TTree *treeR = clustersTree;
dc7ca31d 658
7e08771f 659 AliT0RecPoint frecpoints;
660 AliT0RecPoint * pfrecpoints = &frecpoints;
dc7ca31d 661
662 AliDebug(1,Form("Start FillESD T0"));
663 TBranch *brRec = treeR->GetBranch("T0");
664 if (brRec) {
7e08771f 665 brRec->SetAddress(&pfrecpoints);
dc7ca31d 666 }else{
f16935f7 667 AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d 668 return;
669 }
73df58ab 670
671 brRec->GetEntry(0);
672 Double32_t amp[24], time[24], ampQTC[24], timecorr[24];
673 Double32_t* tcorr;
674 for(Int_t i=0; i<24; i++)
675 amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;
676
ce50812a 677 //1st time
73df58ab 678 Double32_t timeClock[3];
7e08771f 679 Double32_t zPosition = frecpoints.GetVertex();
830172c8 680
ce50812a 681 timeClock[0] = frecpoints.GetT0clock() ;
682 timeClock[1] = frecpoints.Get1stTimeA() + shift;
683 timeClock[2] = frecpoints.Get1stTimeC() - shift;
684 //best time
685 Double32_t timemean[3];
686 timemean[0] = frecpoints.GetMeanTime();
687 timemean[1] = frecpoints.GetBestTimeA() + shift;
688 timemean[2] = frecpoints.GetBestTimeC() - shift;
689
690 for(Int_t i=0; i<3; i++) {
691 fESDTZERO->SetT0TOF(i,timeClock[i]); // interaction time (ns)
692 fESDTZERO->SetT0TOFbest(i,timemean[i]); // interaction time (ns)
693 }
73df58ab 694 for ( Int_t i=0; i<24; i++) {
7e08771f 695 time[i] = frecpoints.GetTime(i); // ps to ns
ce50812a 696 if ( time[i] != 0 && time[i]>-9999) {
7e08771f 697 ampQTC[i] = frecpoints.GetAmp(i);
698 amp[i] = frecpoints.AmpLED(i);
612737bb 699 AliDebug(1,Form("T0: %i time %f ampQTC %f ampLED %f \n", i, time[i], ampQTC[i], amp[i]));
85f61e3b 700 }
73df58ab 701 }
ce50812a 702 fESDTZERO->SetT0time(time); // best TOF on each PMT
703 fESDTZERO->SetT0amplitude(ampQTC); // number of particles(MIPs) on each PMT
7e08771f 704 Int_t trig= frecpoints.GetT0Trig();
705 frecpoints.PrintTriggerSignals( trig);
92f6bfd3 706 // printf(" !!!!! FillESD trigger %i \n",trig);
b0e13b29 707 fESDTZERO->SetT0Trig(trig);
b0e13b29 708 fESDTZERO->SetT0zVertex(zPosition); //vertex Z position
669dc07f 709
7e08771f 710 Double32_t multA=frecpoints.GetMultA();
711 Double32_t multC=frecpoints.GetMultC();
b0e13b29 712 fESDTZERO->SetMultA(multA); // for backward compatubility
713 fESDTZERO->SetMultC(multC); // for backward compatubility
714
715
716 for (Int_t iHit =0; iHit<5; iHit++ ) {
82540e64 717 AliDebug(10,Form("FillESD ::: iHit %i tvdc %f orA %f orC %f\n", iHit,
7e08771f 718 frecpoints.GetTVDC(iHit),
719 frecpoints.GetOrA(iHit),
720 frecpoints.GetOrC(iHit) ));
721 fESDTZERO->SetTVDC(iHit,frecpoints.GetTVDC(iHit));
722 fESDTZERO->SetOrA(iHit,frecpoints.GetOrA(iHit));
723 fESDTZERO->SetOrC(iHit,frecpoints.GetOrC(iHit));
b0e13b29 724
725 for (Int_t i0=0; i0<24; i0++) {
7e08771f 726 // if(frecpoints.GetTimeFull(i0,iHit)>0){
727 // printf("FillESD ::: iHit %i cfd %i time %f \n", iHit, i0, frecpoints.GetTimeFull(i0,iHit));
728 fESDTZERO->SetTimeFull(i0, iHit,frecpoints.GetTimeFull(i0,iHit));
94b4f7e2 729 // }
b0e13b29 730
731 }
732 }
73df58ab 733
ce50812a 734 AliDebug(1,Form("T0: SPDshift %f Vertex %f (T0A+T0C)/2 best %f #ps T0signal %f ps OrA %f ps OrC %f ps T0trig %i\n",shift, zPosition, timemean[0], timeClock[0], timeClock[1], timeClock[2], trig));
b0e13b29 735
ea1a8005 736 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
737 // background flags
738 Bool_t background = BackgroundFlag();
739 fESDTZERO->SetBackgroundFlag(background);
740 Bool_t pileup = PileupFlag();
741 fESDTZERO->SetPileupFlag(pileup);
528890e5 742 for (Int_t i=0; i<5; i++) {
7e08771f 743 fESDTZERO->SetPileupTime(i, frecpoints.GetTVDC(i) ) ;
744 // printf("!!!!!! FillESD :: pileup %i %f %f \n", i,fESDTZERO->GetPileupTime(i), frecpoints.GetTVDC(i));
528890e5 745 }
ea1a8005 746 Bool_t sat = SatelliteFlag();
747 fESDTZERO->SetSatelliteFlag(sat);
73df58ab 748
ea1a8005 749 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
73df58ab 750 if (pESD) {
adf36b9d 751
73df58ab 752 AliESDfriend *fr = (AliESDfriend*)pESD->FindListObject("AliESDfriend");
753 if (fr) {
82540e64 754 AliDebug(10, Form("Writing TZERO friend data to ESD tree"));
73df58ab 755
85f61e3b 756 // if (ncont>2) {
73df58ab 757 tcorr = fESDTZEROfriend->GetT0timeCorr();
758 for ( Int_t i=0; i<24; i++) {
85f61e3b 759 if(i<12 && time[i]>1) timecorr[i] = tcorr[i] - shift/channelWidth;
760 if(i>11 && time[i]>1) timecorr[i] = tcorr[i] + shift/channelWidth;
612737bb 761 if(time[i]>1) AliDebug(10,Form("T0 friend : %i time %f ampQTC %f ampLED %f \n", i, timecorr[i], ampQTC[i], amp[i]));
58bd3a16 762 }
73df58ab 763 fESDTZEROfriend->SetT0timeCorr( timecorr) ;
764 fESDTZEROfriend->SetT0ampLEDminCFD(amp);
765 fESDTZEROfriend->SetT0ampQTC(ampQTC);
85f61e3b 766 fr->SetTZEROfriend(fESDTZEROfriend);
767 // }//
58bd3a16 768 }
b0e13b29 769
770 pESD->SetTZEROData(fESDTZERO);
58bd3a16 771 }
73df58ab 772
dc7ca31d 773} // vertex in 3 sigma
774
ea1a8005 775//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
776 //____________________________________________________________
777
778Bool_t AliT0Reconstructor::PileupFlag() const
779{
780 //
781 Bool_t pileup = false;
528890e5 782 Float_t tvdc[5];
ea1a8005 783 for (Int_t ih=0; ih<5; ih++)
784 {
785 tvdc[ih] = fESDTZERO->GetTVDC(ih);
786
528890e5 787 if( tvdc[0] !=0 && tvdc[0]> -10 && tvdc[0]< 10 )
ea1a8005 788 if(ih>0 && tvdc[ih]>20 ) pileup = true;
528890e5 789 if( tvdc[0] >20 || (tvdc[0] < -20 && tvdc[0] > -9000) ) pileup =true;
790 // if (pileup) printf(" !!!!! pile up %i tvdc %f \n",ih, tvdc[ih]);
ea1a8005 791 }
792
793
794 return pileup;
795
796}
797
798 //____________________________________________________________
799
800Bool_t AliT0Reconstructor::BackgroundFlag() const
801{
802 Bool_t background = false;
803
804 Float_t orA = fESDTZERO->GetOrA(0);
805 Float_t orC = fESDTZERO->GetOrC(0);
806 Float_t tvdc = fESDTZERO->GetTVDC(0);
dc7ca31d 807
ea1a8005 808 if ( (orA > -5 && orA <5) && (orC > -5 && orC <5) && (tvdc < -5 || tvdc > 5)) {
809 background = true;
810 // printf(" orA %f orC %f tvdc %f\n", orA, orC, tvdc);
811 }
812 return background;
813
814
815}
816
817
818 //____________________________________________________________
819
820Bool_t AliT0Reconstructor::SatelliteFlag() const
821{
822
ce50812a 823 Float_t satelliteLow = GetRecoParam() -> GetLowSatelliteThreshold();
824 Float_t satelliteHigh = GetRecoParam() -> GetHighSatelliteThreshold();
ea1a8005 825 Bool_t satellite = false;
826 for (Int_t i0=0; i0<24; i0++) {
827 Float_t timefull = fESDTZERO -> GetTimeFull(i0,0);
ce50812a 828 if( timefull > satelliteLow && timefull < satelliteHigh) satellite=true;
ea1a8005 829 }
830
831 return satellite;
832
833}