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