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