]> 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
Changes for #88334 port to Release code for new T0 reconstruction scheme (Alla)
[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"));
539b9cb9 125
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() ;
ce50812a 147 Double32_t vertex = 9999999, meanVertex = 9999999 ;
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 ) {
196 if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)
d0bcd1fb 197 adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
a7027400 198 else
d0bcd1fb 199 adc[ipmt] = 0;
200
ce50812a 201 time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]), timeCFD->At(ipmt)) ;
202
d0bcd1fb 203 Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
8f620945 204 // time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD->At(ipmt) ) ;
669dc07f 205 AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
ce50812a 206 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
207
644a358e 208 Double_t ampMip = 0;
209 TGraph* ampGraph = (TGraph*)fAmpLED.At(ipmt);
210 if (ampGraph) ampMip = ampGraph->Eval(sl);
211 Double_t qtMip = 0;
212 TGraph* qtGraph = (TGraph*)fQTC.At(ipmt);
213 if (qtGraph) qtMip = qtGraph->Eval(adc[ipmt]);
669dc07f 214 AliDebug(5,Form(" Amlitude in MIPS LED %f , QTC %f in channels %f\n ",ampMip,qtMip, adc[ipmt]));
7e08771f 215 frecpoints.SetTime(ipmt, Float_t(time[ipmt]) );
216 frecpoints.SetAmpLED(ipmt, Float_t( ampMip));
217 frecpoints.SetAmp(ipmt, Float_t(qtMip));
195e1353 218 adcmip[ipmt]=qtMip;
d0bcd1fb 219
dc7ca31d 220 }
221 else {
222 time[ipmt] = 0;
223 adc[ipmt] = 0;
b0ab3f59 224 adcmip[ipmt] = 0;
ce50812a 225
dc7ca31d 226 }
227 }
94c27e4f 228
dc7ca31d 229 for (Int_t ipmt=0; ipmt<12; ipmt++){
ce50812a 230 if(time[ipmt] !=0 && time[ipmt] > -9000
231 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold )
232 {
233 if(time[ipmt]<besttimeC) besttimeC=time[ipmt]; //timeC
234 if(TMath::Abs(time[ipmt])<TMath::Abs(besttimeC_best))
235 besttimeC_best=time[ipmt]; //timeC
dc7ca31d 236 }
dc7ca31d 237 }
ce50812a 238 for ( Int_t ipmt=12; ipmt<24; ipmt++)
239 {
240 if(time[ipmt] != 0 && time[ipmt] > -9000
241 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold)
242 {
243 if(time[ipmt]<besttimeA) besttimeA=time[ipmt];
244 if(TMath::Abs(time[ipmt] ) < TMath::Abs(besttimeA_best))
245 besttimeA_best=time[ipmt]; //timeA
246 }
dc7ca31d 247 }
ce50812a 248
249 if( besttimeA < 999999 && besttimeA!=0) {
250 frecpoints.SetTimeBestA((besttimeA_best * channelWidth - fdZonA/c) - 12500.);
251 frecpoints.SetTime1stA((besttimeA * channelWidth - fdZonA/c) - 12500.);
adf36b9d 252 tr[1]=true;
253 }
ce50812a 254
255 if( besttimeC < 999999 && besttimeC!=0) {
256 frecpoints.SetTimeBestC((besttimeC_best * channelWidth - fdZonC/c) - 12500.);
257 frecpoints.SetTime1stC((besttimeC * channelWidth - fdZonC/c) - 12500.);
adf36b9d 258 tr[2]=true;
259 }
ce50812a 260
261 AliDebug(5,Form("1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f ",
262 besttimeA, besttimeA_best,
263 besttimeC, besttimeC_best) );
264
adf36b9d 265 if(besttimeA <999999 && besttimeC < 999999 ){
9b83615d 266 // timeDiff = (besttimeC - besttimeA)*channelWidth;
267 timeDiff = (besttimeA - besttimeC)*channelWidth;
ce50812a 268 meanTime = channelWidth * (besttimeA_best + besttimeC_best)/2. ;
269 timeclock = channelWidth * (besttimeA + besttimeC)/2. ;
adf36b9d 270 vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
271 tr[0]=true;
272 }
7e08771f 273 frecpoints.SetVertex(vertex);
ce50812a 274 frecpoints.SetMeanTime(meanTime - 12500.);
275 frecpoints.SetT0clock(timeclock - 12500.);
7e08771f 276 frecpoints.SetT0Trig(tr);
ce50812a 277
278 AliDebug(5,Form("fRecPoints::: 1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f shiftA %f shiftC %f ",
279 frecpoints.Get1stTimeA(), frecpoints.GetBestTimeA(),
280 frecpoints.Get1stTimeC(), frecpoints.GetBestTimeC(),
281 fTimeMeanShift[1],fTimeMeanShift[2] ) );
282
669dc07f 283 AliDebug(5,Form("T0 triggers %d %d %d %d %d",tr[0],tr[1],tr[2],tr[3],tr[4]));
ce50812a 284
adf36b9d 285 //online mean
7e08771f 286 frecpoints.SetOnlineMean(Int_t(onlineMean));
b5a9f753 287 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 288
830172c8 289 clustersTree->Fill();
ce50812a 290
bd375212 291 delete timeCFD;
292 delete timeLED;
293 delete chargeQT0;
294 delete chargeQT1;
dc7ca31d 295}
296
297
c41ceaac 298//_______________________________________________________________________
299
300void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
301{
94c27e4f 302 // T0 raw ->
539b9cb9 303 //
b0e13b29 304
305 Float_t meanOrA = fTime0vertex[0] + 587;
306 Float_t meanOrC = fTime0vertex[0] + 678;
307 Float_t meanTVDC = fTime0vertex[0] + 2564;
82540e64 308 Int_t timeDelayCFD[24];
309
d3e04608 310
612737bb 311 Int_t badpmt[24];
38cbfa7c 312 //Bad channel
58432641 313 for (Int_t i=0; i<24; i++) {
314 badpmt[i] = GetRecoParam() -> GetBadChannels(i);
82540e64 315 timeDelayCFD[i] = Int_t (fParam->GetTimeDelayCFD(i));
58432641 316 }
ce50812a 317 Int_t equalize = GetRecoParam() -> GetEq();
318 printf( "AliT0Reconstructor::Reconstruct::: RecoParam %i \n",equalize);
319 fCalib->SetEq(equalize);
d3e04608 320 Int_t low[500], high[500];
94b4f7e2 321 Float_t timefull=-9999;;
322 Float_t tvdc = -9999; Float_t ora = -9999; Float_t orc = -9999;
58bd3a16 323
e8ed1cd0 324 Int_t allData[110][5];
2e6a5ee0 325
e8ed1cd0 326 Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
ce50812a 327 Float_t time2zero[24];
612737bb 328 Double32_t timeDiff, meanTime, timeclock;
329 timeDiff = meanTime = timeclock = 9999999;
adf36b9d 330 Float_t c = 29.9792458; // cm/ns
331 Double32_t vertex = 9999999;
776de217 332 Int_t onlineMean=0;
9480f05f 333 Float_t meanVertex = 0;
612737bb 334 for (Int_t i0=0; i0<24; i0++) {
335 low[i0] = Int_t(fTime0vertex[i0]) - 200;
336 high[i0] = Int_t(fTime0vertex[i0]) + 200;
ce50812a 337 time2zero[i0] = 99999;
612737bb 338 }
339
669dc07f 340 for (Int_t i0=0; i0<110; i0++)
b0e13b29 341 for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
612737bb 342
343 Float_t lowAmpThreshold = GetRecoParam()->GetAmpLowThreshold();
344 Float_t highAmpThreshold = GetRecoParam()->GetAmpHighThreshold();
ce50812a 345 printf( "AliT0Reconstructor::Reconstruct::: RecoParam amplitude %f %f \n",lowAmpThreshold, highAmpThreshold);
346
347 Double32_t besttimeA=9999999; Double32_t besttimeA_best=9999999;
348 Double32_t besttimeC=9999999; Double32_t besttimeC_best=9999999;
349
350 Float_t channelWidth = fParam->GetChannelWidth() ;
b0e13b29 351
7e08771f 352 AliT0RecPoint frecpoints;
353 AliT0RecPoint * pfrecpoints = &frecpoints;
bce12dc5 354
7e08771f 355 recTree->Branch( "T0", "AliT0RecPoint" ,&pfrecpoints);
2e6a5ee0 356
bce12dc5 357 AliDebug(10," before read data ");
358 AliT0RawReader myrawreader(rawReader);
776de217 359
360 UInt_t type =rawReader->GetType();
361
bce12dc5 362 if (!myrawreader.Next())
363 AliDebug(1,Form(" no raw data found!!"));
364 else
365 {
38cbfa7c 366 for (Int_t i=0; i<24; i++)
367 {
368 timeCFD[i]=0; timeLED[i]=0; chargeQT0[i]=0; chargeQT1[i]=0;
369 }
370 Int_t fBCID=Int_t (rawReader->GetBCID());
8f620945 371 Int_t trmbunch= myrawreader.GetTRMBunchID();
86fd0587 372 AliDebug(10,Form(" CDH BC ID %i, TRM BC ID %i \n", fBCID, trmbunch ));
373
38cbfa7c 374 if(type == 7 ) { //only physics
669dc07f 375 for (Int_t i=0; i<107; i++) {
bce12dc5 376 for (Int_t iHit=0; iHit<5; iHit++)
377 {
378 allData[i][iHit] = myrawreader.GetData(i,iHit);
379 }
8f620945 380 }
8f620945 381
85f61e3b 382 for (Int_t in=0; in<12; in++)
383 {
384 for (Int_t iHit=0; iHit<5; iHit++)
38cbfa7c 385 {
612737bb 386 if(allData[in+1][iHit] > low[in] &&
387 allData[in+1][iHit] < high[in])
38cbfa7c 388 {
85f61e3b 389 timeCFD[in] = allData[in+1][iHit] ;
390 break;
38cbfa7c 391 }
85f61e3b 392 }
393 for (Int_t iHit=0; iHit<5; iHit++)
394 {
612737bb 395 if(allData[in+1+56][iHit] > low[in] &&
396 allData[in+1+56][iHit] < high[in])
38cbfa7c 397 {
85f61e3b 398 timeCFD[in+12] = allData[in+56+1][iHit] ;
399 break;
38cbfa7c 400 }
38cbfa7c 401 }
612737bb 402 timeLED[in+12] = allData[in+68+1][0] ;
403 timeLED[in] = allData[in+12+1][0] ;
ce50812a 404 AliDebug(50, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
85f61e3b 405 in, timeCFD[in],timeCFD[in+12],timeLED[in],
612737bb 406 timeLED[in+12]));
38cbfa7c 407
8f620945 408 }
409
8f620945 410
85f61e3b 411 for (Int_t in=0; in<12; in++)
412 {
413 chargeQT0[in]=allData[2*in+25][0];
414 chargeQT1[in]=allData[2*in+26][0];
612737bb 415 AliDebug(25, Form(" readed Raw %i %i %i",
85f61e3b 416 in, chargeQT0[in],chargeQT1[in]));
417 }
418 for (Int_t in=12; in<24; in++)
419 {
420 chargeQT0[in]=allData[2*in+57][0];
421 chargeQT1[in]=allData[2*in+58][0];
612737bb 422 AliDebug(25, Form(" readed Raw %i %i %i",
85f61e3b 423 in, chargeQT0[in],chargeQT1[in]));
85f61e3b 424 }
425
8f226345 426 onlineMean = allData[49][0];
427
195e1353 428 Double32_t time[24], adc[24], adcmip[24], noncalibtime[24];
8f620945 429 for (Int_t ipmt=0; ipmt<24; ipmt++) {
36bfca7d 430 if(timeCFD[ipmt] > 0 /* && badpmt[ipmt]==0*/ ){
bce12dc5 431 //for simulated data
d0bcd1fb 432 //for physics data
73df58ab 433 if(( chargeQT0[ipmt] - chargeQT1[ipmt])>0) {
541b42c4 434 adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
73df58ab 435 }
d0bcd1fb 436 else
437 adc[ipmt] = 0;
612737bb 438 // time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
bce12dc5 439
612737bb 440 time[ipmt] = fCalib-> WalkCorrection(Int_t (fTime0vertex[ipmt]), ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
d0bcd1fb 441 Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
8f620945 442 // time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
73df58ab 443 AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
d0bcd1fb 444 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
644a358e 445 Double_t ampMip = 0;
446 TGraph * ampGraph = (TGraph*)fAmpLED.At(ipmt);
447 if (ampGraph) ampMip = ampGraph->Eval(sl);
448 Double_t qtMip = 0;
449 TGraph * qtGraph = (TGraph*)fQTC.At(ipmt);
450 if (qtGraph) qtMip = qtGraph->Eval(adc[ipmt]);
345f03db 451 AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %f\n ",ampMip,qtMip, adc[ipmt]));
ce50812a 452 if( equalize ==0 )
7e08771f 453 frecpoints.SetTime(ipmt, Float_t(time[ipmt]) );
ce50812a 454 else
455 frecpoints.SetTime(ipmt, Float_t(time[ipmt] + fTime0vertex[ipmt]) );
456
457 frecpoints.SetAmp(ipmt, Double32_t( qtMip));
458 adcmip[ipmt]=qtMip;
459 frecpoints.SetAmpLED(ipmt, Double32_t(ampMip));
460 noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
2e6a5ee0 461 }
462 else {
ce50812a 463 time[ipmt] = -9999;
2e6a5ee0 464 adc[ipmt] = 0;
b0ab3f59 465 adcmip[ipmt] = 0;
ce50812a 466 noncalibtime[ipmt] = -9999;
2e6a5ee0 467 }
468 }
ce50812a 469 fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;
470
2e6a5ee0 471 for (Int_t ipmt=0; ipmt<12; ipmt++){
ce50812a 472 if(time[ipmt] !=0 && time[ipmt] > -9000
473 /*&& badpmt[ipmt]==0 */
474 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold )
73df58ab 475 {
ce50812a 476 if(time[ipmt]<besttimeC) besttimeC=time[ipmt]; //timeC
477 if(TMath::Abs(time[ipmt])<TMath::Abs(besttimeC_best))
478 besttimeC_best=time[ipmt]; //timeC
2e6a5ee0 479 }
2e6a5ee0 480 }
73df58ab 481 for ( Int_t ipmt=12; ipmt<24; ipmt++)
482 {
ce50812a 483 if(time[ipmt] != 0 && time[ipmt] > -9000
484 /* && badpmt[ipmt]==0*/
485 && adcmip[ipmt]>lowAmpThreshold && adcmip[ipmt]<highAmpThreshold)
73df58ab 486 {
ce50812a 487 if(time[ipmt]<besttimeA) besttimeA=time[ipmt];
488 if(TMath::Abs(time[ipmt] ) < TMath::Abs(besttimeA_best))
489 besttimeA_best=time[ipmt]; //timeA
73df58ab 490 }
2e6a5ee0 491 }
ce50812a 492
493 if(besttimeA < 999999 && besttimeA!=0 ) {
494 if( equalize ==0 )
495 frecpoints.SetTime1stA((besttimeA * channelWidth)- 1000.*fLatencyHPTDC + 1000.*fLatencyL1A - 1000.*fGRPdelays - fTimeMeanShift[1] );
496 else
497 {
498 frecpoints.SetTimeBestA((besttimeA_best * channelWidth ));
499 frecpoints.SetTime1stA((besttimeA * channelWidth - fTimeMeanShift[1]));
500 }
501 }
502 if( besttimeC < 999999 && besttimeC!=0) {
503 if( equalize ==0 )
504 frecpoints.SetTime1stC((besttimeC * channelWidth)- 1000.*fLatencyHPTDC +1000.*fLatencyL1C - 1000.*fGRPdelays - fTimeMeanShift[2]);
505 else
506 {
507 frecpoints.SetTimeBestC((besttimeC_best * channelWidth ));
508 frecpoints.SetTime1stC((besttimeC * channelWidth - fTimeMeanShift[2]));
509 }
510 }
511 AliDebug(5,Form("1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f ",
512 besttimeA, besttimeA_best,
513 besttimeC, besttimeC_best) );
514 AliDebug(5,Form("fRecPoints::: 1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f shiftA %f shiftC %f ",
515 frecpoints.Get1stTimeA(), frecpoints.GetBestTimeA(),
516 frecpoints.Get1stTimeC(), frecpoints.GetBestTimeC(),
517 fTimeMeanShift[1],fTimeMeanShift[2] ) );
518 if( besttimeC < 999999 && besttimeA < 999999) {
519 timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - fTimeMeanShift[0];
520 meanTime = channelWidth * Float_t(besttimeA_best + besttimeC_best )/2.;
612737bb 521 timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth ;
9480f05f 522 vertex = meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2;
ce50812a 523 }
524
776de217 525 } //if phys event
ce50812a 526 AliDebug(1,Form(" timeDiff %f #channel, meanTime %f #ps, TOFmean%f vertex %f cm meanVertex %f \n",timeDiff, meanTime,timeclock, vertex,meanVertex));
527 frecpoints.SetT0clock(timeclock);
528 frecpoints.SetVertex(vertex);
529 frecpoints.SetMeanTime(meanTime);
530 frecpoints.SetOnlineMean(Int_t(onlineMean));
adf36b9d 531
ce50812a 532 // Set triggers
adf36b9d 533 Bool_t tr[5];
92f6bfd3 534 Int_t trchan[5] = {50,51,52,55,56};
535 Float_t lowtr[5] = {meanTVDC-700, meanOrA-700, meanOrC-700, meanOrC-1000, meanOrC-1000 };
536 Float_t hightr[5] = {meanTVDC+700, meanOrA+700, meanOrC+700, meanOrC+1000, meanOrC+1000};
ce50812a 537
82540e64 538 for (Int_t i=0; i<5; i++) tr[i] = false;
adf36b9d 539 for (Int_t itr=0; itr<5; itr++) {
b0e13b29 540 for (Int_t iHit=0; iHit<1; iHit++)
38cbfa7c 541 {
542 Int_t trr=trchan[itr];
92f6bfd3 543 if( allData[trr][iHit] > lowtr[itr] && allData[trr][iHit] < hightr[itr]) tr[itr]=true;
ce50812a 544
545 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 546 }
adf36b9d 547 }
7e08771f 548 frecpoints.SetT0Trig(tr);
b0e13b29 549
ce50812a 550 // all times with amplitude correction
551 Float_t timecent;
b0e13b29 552 for (Int_t iHit=0; iHit<5; iHit++)
553 {
ce50812a 554 timefull = timecent = -9999;
82540e64 555 tvdc = ora = orc = -9999;
b0e13b29 556 if(allData[50][iHit]>0)
557 tvdc = (Float_t(allData[50][iHit]) - meanTVDC) * channelWidth* 0.001;
558 if(allData[51][iHit]>0)
559 ora = (Float_t(allData[51][iHit]) - meanOrA) * channelWidth* 0.001;
560
561 if(allData[52][iHit]>0)
562 orc = (Float_t(allData[52][iHit]) - meanOrC) * channelWidth* 0.001;
563
7e08771f 564 frecpoints.SetOrC( iHit, orc);
565 frecpoints.SetOrA( iHit, ora);
566 frecpoints.SetTVDC( iHit, tvdc);
b0e13b29 567 for (Int_t i0=0; i0<12; i0++) {
ce50812a 568 if (equalize ==0 )
569 timecent = fTime0vertex[i0] + timeDelayCFD[i0];
570 else
571 timecent = fTime0vertex[i0];
b0e13b29 572 timefull = -9999;
573 if(allData[i0+1][iHit]>1)
ce50812a 574 timefull = (Float_t(allData[i0+1][iHit]) - timecent)* channelWidth* 0.001;
7e08771f 575 frecpoints.SetTimeFull(i0, iHit,timefull) ;
82540e64 576 // 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 577
578 }
579
580 for (Int_t i0=12; i0<24; i0++) {
581 timefull = -9999;
ce50812a 582 if (equalize ==0 )
583 timecent = fTime0vertex[i0] + timeDelayCFD[i0];
584 else
585 timecent = fTime0vertex[i0];
94b4f7e2 586 if(allData[i0+45][iHit]>1) {
ce50812a 587 timefull = (Float_t(allData[i0+45][iHit]) - timecent)* channelWidth* 0.001;
94b4f7e2 588 }
82540e64 589 // 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 590 frecpoints.SetTimeFull(i0, iHit, timefull) ;
b0e13b29 591 }
592 }
593
594
669dc07f 595 //Set MPD
596 if(allData[53][0]>0 && allData[54][0])
7e08771f 597 frecpoints.SetMultA(allData[53][0]-allData[54][0]);
b0e13b29 598 if(allData[105][0]>0 && allData[106][0])
7e08771f 599 frecpoints.SetMultC(allData[105][0]-allData[106][0]);
b0e13b29 600
601
669dc07f 602 } // if (else )raw data
58bd3a16 603 recTree->Fill();
58bd3a16 604}
adf36b9d 605
606
ce50812a 607//____________________________________________________________
adf36b9d 608
609 void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
610 {
dc7ca31d 611
612 /***************************************************
613 Resonstruct digits to vertex position
614 ****************************************************/
615
dc7ca31d 616 AliDebug(1,Form("Start FillESD T0"));
b0ab3f59 617 if(!pESD) {
618 AliError("No ESD Event");
619 return;
620 }
4cbe597e 621 pESD ->SetT0spread(fTimeSigmaShift);
622
36cde487 623
58bd3a16 624 Float_t channelWidth = fParam->GetChannelWidth() ;
f7c2c2fc 625 Float_t c = 0.0299792458; // cm/ps
adf36b9d 626 Float_t currentVertex=0, shift=0;
291f31a1 627 Int_t ncont=-1;
adf36b9d 628 const AliESDVertex* vertex = pESD->GetPrimaryVertex();
629 if (!vertex) vertex = pESD->GetPrimaryVertexSPD();
630 if (!vertex) vertex = pESD->GetPrimaryVertexTPC();
631 if (!vertex) vertex = pESD->GetVertex();
632
633 if (vertex) {
634 AliDebug(2, Form("Got %s (%s) from ESD: %f",
635 vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
636 currentVertex = vertex->GetZ();
637
638 ncont = vertex->GetNContributors();
291f31a1 639 if(ncont>0 ) {
adf36b9d 640 shift = currentVertex/c;
adf36b9d 641 }
642 }
d76c31f4 643 TTree *treeR = clustersTree;
dc7ca31d 644
7e08771f 645 AliT0RecPoint frecpoints;
646 AliT0RecPoint * pfrecpoints = &frecpoints;
dc7ca31d 647
648 AliDebug(1,Form("Start FillESD T0"));
649 TBranch *brRec = treeR->GetBranch("T0");
650 if (brRec) {
7e08771f 651 brRec->SetAddress(&pfrecpoints);
dc7ca31d 652 }else{
f16935f7 653 AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d 654 return;
655 }
73df58ab 656
657 brRec->GetEntry(0);
658 Double32_t amp[24], time[24], ampQTC[24], timecorr[24];
659 Double32_t* tcorr;
660 for(Int_t i=0; i<24; i++)
661 amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;
662
ce50812a 663 //1st time
73df58ab 664 Double32_t timeClock[3];
7e08771f 665 Double32_t zPosition = frecpoints.GetVertex();
830172c8 666
ce50812a 667 timeClock[0] = frecpoints.GetT0clock() ;
668 timeClock[1] = frecpoints.Get1stTimeA() + shift;
669 timeClock[2] = frecpoints.Get1stTimeC() - shift;
670 //best time
671 Double32_t timemean[3];
672 timemean[0] = frecpoints.GetMeanTime();
673 timemean[1] = frecpoints.GetBestTimeA() + shift;
674 timemean[2] = frecpoints.GetBestTimeC() - shift;
675
676 for(Int_t i=0; i<3; i++) {
677 fESDTZERO->SetT0TOF(i,timeClock[i]); // interaction time (ns)
678 fESDTZERO->SetT0TOFbest(i,timemean[i]); // interaction time (ns)
679 }
73df58ab 680 for ( Int_t i=0; i<24; i++) {
7e08771f 681 time[i] = frecpoints.GetTime(i); // ps to ns
ce50812a 682 if ( time[i] != 0 && time[i]>-9999) {
7e08771f 683 ampQTC[i] = frecpoints.GetAmp(i);
684 amp[i] = frecpoints.AmpLED(i);
612737bb 685 AliDebug(1,Form("T0: %i time %f ampQTC %f ampLED %f \n", i, time[i], ampQTC[i], amp[i]));
85f61e3b 686 }
73df58ab 687 }
ce50812a 688 fESDTZERO->SetT0time(time); // best TOF on each PMT
689 fESDTZERO->SetT0amplitude(ampQTC); // number of particles(MIPs) on each PMT
7e08771f 690 Int_t trig= frecpoints.GetT0Trig();
691 frecpoints.PrintTriggerSignals( trig);
92f6bfd3 692 // printf(" !!!!! FillESD trigger %i \n",trig);
b0e13b29 693 fESDTZERO->SetT0Trig(trig);
b0e13b29 694 fESDTZERO->SetT0zVertex(zPosition); //vertex Z position
669dc07f 695
7e08771f 696 Double32_t multA=frecpoints.GetMultA();
697 Double32_t multC=frecpoints.GetMultC();
b0e13b29 698 fESDTZERO->SetMultA(multA); // for backward compatubility
699 fESDTZERO->SetMultC(multC); // for backward compatubility
700
701
702 for (Int_t iHit =0; iHit<5; iHit++ ) {
82540e64 703 AliDebug(10,Form("FillESD ::: iHit %i tvdc %f orA %f orC %f\n", iHit,
7e08771f 704 frecpoints.GetTVDC(iHit),
705 frecpoints.GetOrA(iHit),
706 frecpoints.GetOrC(iHit) ));
707 fESDTZERO->SetTVDC(iHit,frecpoints.GetTVDC(iHit));
708 fESDTZERO->SetOrA(iHit,frecpoints.GetOrA(iHit));
709 fESDTZERO->SetOrC(iHit,frecpoints.GetOrC(iHit));
b0e13b29 710
711 for (Int_t i0=0; i0<24; i0++) {
7e08771f 712 // if(frecpoints.GetTimeFull(i0,iHit)>0){
713 // printf("FillESD ::: iHit %i cfd %i time %f \n", iHit, i0, frecpoints.GetTimeFull(i0,iHit));
714 fESDTZERO->SetTimeFull(i0, iHit,frecpoints.GetTimeFull(i0,iHit));
94b4f7e2 715 // }
b0e13b29 716
717 }
718 }
73df58ab 719
ce50812a 720 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 721
ea1a8005 722 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
723 // background flags
724 Bool_t background = BackgroundFlag();
725 fESDTZERO->SetBackgroundFlag(background);
726 Bool_t pileup = PileupFlag();
727 fESDTZERO->SetPileupFlag(pileup);
528890e5 728 for (Int_t i=0; i<5; i++) {
7e08771f 729 fESDTZERO->SetPileupTime(i, frecpoints.GetTVDC(i) ) ;
730 // printf("!!!!!! FillESD :: pileup %i %f %f \n", i,fESDTZERO->GetPileupTime(i), frecpoints.GetTVDC(i));
528890e5 731 }
ea1a8005 732 Bool_t sat = SatelliteFlag();
733 fESDTZERO->SetSatelliteFlag(sat);
73df58ab 734
ea1a8005 735 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
73df58ab 736 if (pESD) {
adf36b9d 737
73df58ab 738 AliESDfriend *fr = (AliESDfriend*)pESD->FindListObject("AliESDfriend");
739 if (fr) {
82540e64 740 AliDebug(10, Form("Writing TZERO friend data to ESD tree"));
73df58ab 741
85f61e3b 742 // if (ncont>2) {
73df58ab 743 tcorr = fESDTZEROfriend->GetT0timeCorr();
744 for ( Int_t i=0; i<24; i++) {
85f61e3b 745 if(i<12 && time[i]>1) timecorr[i] = tcorr[i] - shift/channelWidth;
746 if(i>11 && time[i]>1) timecorr[i] = tcorr[i] + shift/channelWidth;
612737bb 747 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 748 }
73df58ab 749 fESDTZEROfriend->SetT0timeCorr( timecorr) ;
750 fESDTZEROfriend->SetT0ampLEDminCFD(amp);
751 fESDTZEROfriend->SetT0ampQTC(ampQTC);
85f61e3b 752 fr->SetTZEROfriend(fESDTZEROfriend);
753 // }//
58bd3a16 754 }
b0e13b29 755
756 pESD->SetTZEROData(fESDTZERO);
58bd3a16 757 }
73df58ab 758
dc7ca31d 759} // vertex in 3 sigma
760
ea1a8005 761//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
762 //____________________________________________________________
763
764Bool_t AliT0Reconstructor::PileupFlag() const
765{
766 //
767 Bool_t pileup = false;
528890e5 768 Float_t tvdc[5];
ea1a8005 769 for (Int_t ih=0; ih<5; ih++)
770 {
771 tvdc[ih] = fESDTZERO->GetTVDC(ih);
772
528890e5 773 if( tvdc[0] !=0 && tvdc[0]> -10 && tvdc[0]< 10 )
ea1a8005 774 if(ih>0 && tvdc[ih]>20 ) pileup = true;
528890e5 775 if( tvdc[0] >20 || (tvdc[0] < -20 && tvdc[0] > -9000) ) pileup =true;
776 // if (pileup) printf(" !!!!! pile up %i tvdc %f \n",ih, tvdc[ih]);
ea1a8005 777 }
778
779
780 return pileup;
781
782}
783
784 //____________________________________________________________
785
786Bool_t AliT0Reconstructor::BackgroundFlag() const
787{
788 Bool_t background = false;
789
790 Float_t orA = fESDTZERO->GetOrA(0);
791 Float_t orC = fESDTZERO->GetOrC(0);
792 Float_t tvdc = fESDTZERO->GetTVDC(0);
dc7ca31d 793
ea1a8005 794 if ( (orA > -5 && orA <5) && (orC > -5 && orC <5) && (tvdc < -5 || tvdc > 5)) {
795 background = true;
796 // printf(" orA %f orC %f tvdc %f\n", orA, orC, tvdc);
797 }
798 return background;
799
800
801}
802
803
804 //____________________________________________________________
805
806Bool_t AliT0Reconstructor::SatelliteFlag() const
807{
808
ce50812a 809 Float_t satelliteLow = GetRecoParam() -> GetLowSatelliteThreshold();
810 Float_t satelliteHigh = GetRecoParam() -> GetHighSatelliteThreshold();
ea1a8005 811 Bool_t satellite = false;
812 for (Int_t i0=0; i0<24; i0++) {
813 Float_t timefull = fESDTZERO -> GetTimeFull(i0,0);
ce50812a 814 if( timefull > satelliteLow && timefull < satelliteHigh) satellite=true;
ea1a8005 815 }
816
817 return satellite;
818
819}