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