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