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