removed huge printout
[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();
993257ce 87
85f61e3b 88 fGRPdelays = l1Delay - phase->GetMeanPhase();
993257ce 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();
993257ce 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);
612737bb 110 }
58bd3a16 111 fLatencyL1 = fParam->GetLatencyL1();
612737bb 112 fLatencyL1A = fParam->GetLatencyL1A();
58bd3a16 113 fLatencyL1C = fParam->GetLatencyL1C();
114 fLatencyHPTDC = fParam->GetLatencyHPTDC();
9d026202 115 AliDebug(2,Form(" LatencyL1 %f latencyL1A %f latencyL1C %f latencyHPTDC %f \n",fLatencyL1, fLatencyL1A, fLatencyL1C, fLatencyHPTDC));
612737bb 116
117 for (Int_t i=0; i<24; i++) {
993257ce 118 if( fTime0vertex[i] < 500 || fTime0vertex[i] > 60000) fTime0vertex[i] =( 1000.*fLatencyHPTDC - 1000.*fLatencyL1 + 1000.*fGRPdelays)/24.4;
119 AliDebug(2,Form("OCDB mean CFD time %i %f \n",i, fTime0vertex[i]));
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"));
993257ce 124
12e9daf9 125 fCalib = new AliT0Calibrator();
73df58ab 126 fESDTZEROfriend = new AliESDTZEROfriend();
b0e13b29 127 fESDTZERO = new AliESDTZERO();
993257ce 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);
993257ce 160 AliDebug(2, Form("Reconstruct(TTree*digitsTree shiftA %f shiftC %f shiftAC %f \n",shiftA, shiftC, shiftAC));
161
ce50812a 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++){
993257ce 237 if(time[ipmt] !=0 && time[ipmt] > 0
ce50812a 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 {
993257ce 247 if(time[ipmt] != 0 && time[ipmt] > 0
ce50812a 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 //
993257ce 311
b0e13b29 312 Float_t meanOrA = fTime0vertex[0] + 587;
313 Float_t meanOrC = fTime0vertex[0] + 678;
314 Float_t meanTVDC = fTime0vertex[0] + 2564;
993257ce 315 Float_t meanQT1 = fTime0vertex[0] + 2564;
316 Float_t meanQT0 = fTime0vertex[0] + 3564;
317
82540e64 318 Int_t timeDelayCFD[24];
5773ee77 319 Int_t corridor = GetRecoParam() -> GetCorridor();
13e2fbbd 320 // printf("!!!! corridor %i \n",corridor);
612737bb 321 Int_t badpmt[24];
38cbfa7c 322 //Bad channel
58432641 323 for (Int_t i=0; i<24; i++) {
324 badpmt[i] = GetRecoParam() -> GetBadChannels(i);
82540e64 325 timeDelayCFD[i] = Int_t (fParam->GetTimeDelayCFD(i));
58432641 326 }
ce50812a 327 Int_t equalize = GetRecoParam() -> GetEq();
13e2fbbd 328 printf( "AliT0Reconstructor::Reconstruct::: RecoParam %i \n",equalize);
ce50812a 329 fCalib->SetEq(equalize);
d3e04608 330 Int_t low[500], high[500];
13e2fbbd 331 Float_t timefull=-99999;;
332 Float_t tvdc = -99999; Float_t ora = -99999; Float_t orc = -99999;
993257ce 333
e8ed1cd0 334 Int_t allData[110][5];
2e6a5ee0 335
e8ed1cd0 336 Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
ce50812a 337 Float_t time2zero[24];
612737bb 338 Double32_t timeDiff, meanTime, timeclock;
339 timeDiff = meanTime = timeclock = 9999999;
adf36b9d 340 Float_t c = 29.9792458; // cm/ns
341 Double32_t vertex = 9999999;
776de217 342 Int_t onlineMean=0;
9480f05f 343 Float_t meanVertex = 0;
b7fd8e80 344 Int_t pedestal[24];
612737bb 345 for (Int_t i0=0; i0<24; i0++) {
5773ee77 346 low[i0] = Int_t(fTime0vertex[i0]) - corridor;
347 high[i0] = Int_t(fTime0vertex[i0]) + corridor;
348 time2zero[i0] = 99999;
349 pedestal[i0]=Int_t (GetRecoParam()->GetLow(100+i0) );
612737bb 350 }
351
669dc07f 352 for (Int_t i0=0; i0<110; i0++)
b0e13b29 353 for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
612737bb 354
355 Float_t lowAmpThreshold = GetRecoParam()->GetAmpLowThreshold();
356 Float_t highAmpThreshold = GetRecoParam()->GetAmpHighThreshold();
fd3fb735 357
ce50812a 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
85f61e3b 425 for (Int_t in=0; in<12; in++)
426 {
993257ce 427 for (Int_t iHit=0; iHit<5; iHit++)
428 {
429 if (allData[2*in+26][iHit] > fTime0vertex[0]+2000 &&
430 allData[2*in+26][iHit] <fTime0vertex[0]+3000 &&
431 allData[2*in+25][iHit] > fTime0vertex[0]+3000)
794ab872 432 {
794ab872 433 chargeQT0[in]=allData[2*in+25][0];
434 chargeQT1[in]=allData[2*in+26][0];
435 AliDebug(25, Form(" readed Raw %i %i %i",
436 in, chargeQT0[in],chargeQT1[in]));
993257ce 437 break;
794ab872 438 }
993257ce 439 }
85f61e3b 440 }
441 for (Int_t in=12; in<24; in++)
442 {
993257ce 443 for (Int_t iHit=0; iHit<5; iHit++)
444 {
445 if (allData[2*in+58][iHit] > fTime0vertex[0]+2000 &&
446 allData[2*in+58][iHit] <fTime0vertex[0]+3000 &&
447 allData[2*in+57][iHit] > fTime0vertex[0]+3000)
794ab872 448 {
449 chargeQT0[in]=allData[2*in+57][0];
450 chargeQT1[in]=allData[2*in+58][0];
451 AliDebug(25, Form(" readed Raw %i %i %i",
452 in, chargeQT0[in],chargeQT1[in]));
993257ce 453 break;
454 }
794ab872 455 }
85f61e3b 456 }
457
8f226345 458 onlineMean = allData[49][0];
459
195e1353 460 Double32_t time[24], adc[24], adcmip[24], noncalibtime[24];
8f620945 461 for (Int_t ipmt=0; ipmt<24; ipmt++) {
5773ee77 462 if(timeCFD[ipmt] > 0 && (chargeQT0[ipmt] - chargeQT1[ipmt])> 0 ){
463 //for simulated data
d0bcd1fb 464 //for physics data
993257ce 465 adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
c8d939c8 466 Int_t refAmp = Int_t (fTime0vertex[ipmt]);
467 time[ipmt] = fCalib-> WalkCorrection( refAmp, ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
d0bcd1fb 468 Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
8f620945 469 // time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
73df58ab 470 AliDebug(5,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
d0bcd1fb 471 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
644a358e 472 Double_t ampMip = 0;
473 TGraph * ampGraph = (TGraph*)fAmpLED.At(ipmt);
474 if (ampGraph) ampMip = ampGraph->Eval(sl);
475 Double_t qtMip = 0;
476 TGraph * qtGraph = (TGraph*)fQTC.At(ipmt);
477 if (qtGraph) qtMip = qtGraph->Eval(adc[ipmt]);
345f03db 478 AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %f\n ",ampMip,qtMip, adc[ipmt]));
13e2fbbd 479 if( equalize ==0 )
480 frecpoints.SetTime(ipmt, Float_t(time[ipmt]) );
481 else
482 frecpoints.SetTime(ipmt, Float_t(time[ipmt] + fTime0vertex[ipmt]) );
483 // frecpoints.SetTime(ipmt, Float_t(time[ipmt] ) );
ce50812a 484 frecpoints.SetAmp(ipmt, Double32_t( qtMip));
485 adcmip[ipmt]=qtMip;
5773ee77 486 frecpoints.SetAmpLED(ipmt, Double32_t(ampMip));
ce50812a 487 noncalibtime[ipmt]= Double32_t (timeCFD[ipmt]);
13e2fbbd 488 }
489 else {
ce50812a 490 time[ipmt] = -9999;
2e6a5ee0 491 adc[ipmt] = 0;
b0ab3f59 492 adcmip[ipmt] = 0;
ce50812a 493 noncalibtime[ipmt] = -9999;
13e2fbbd 494 }
2e6a5ee0 495 }
ce50812a 496 fESDTZEROfriend->SetT0timeCorr(noncalibtime) ;
497
2e6a5ee0 498 for (Int_t ipmt=0; ipmt<12; ipmt++){
ce50812a 499 if(time[ipmt] !=0 && time[ipmt] > -9000
500 /*&& badpmt[ipmt]==0 */
13e2fbbd 501 && adcmip[ipmt]>lowAmpThreshold )
73df58ab 502 {
ce50812a 503 if(time[ipmt]<besttimeC) besttimeC=time[ipmt]; //timeC
504 if(TMath::Abs(time[ipmt])<TMath::Abs(besttimeC_best))
505 besttimeC_best=time[ipmt]; //timeC
2e6a5ee0 506 }
2e6a5ee0 507 }
73df58ab 508 for ( Int_t ipmt=12; ipmt<24; ipmt++)
509 {
ce50812a 510 if(time[ipmt] != 0 && time[ipmt] > -9000
511 /* && badpmt[ipmt]==0*/
13e2fbbd 512 && adcmip[ipmt]>lowAmpThreshold )
73df58ab 513 {
ce50812a 514 if(time[ipmt]<besttimeA) besttimeA=time[ipmt];
515 if(TMath::Abs(time[ipmt] ) < TMath::Abs(besttimeA_best))
516 besttimeA_best=time[ipmt]; //timeA
73df58ab 517 }
2e6a5ee0 518 }
ce50812a 519
520 if(besttimeA < 999999 && besttimeA!=0 ) {
521 if( equalize ==0 )
522 frecpoints.SetTime1stA((besttimeA * channelWidth)- 1000.*fLatencyHPTDC + 1000.*fLatencyL1A - 1000.*fGRPdelays - fTimeMeanShift[1] );
523 else
524 {
525 frecpoints.SetTimeBestA((besttimeA_best * channelWidth ));
526 frecpoints.SetTime1stA((besttimeA * channelWidth - fTimeMeanShift[1]));
527 }
528 }
529 if( besttimeC < 999999 && besttimeC!=0) {
530 if( equalize ==0 )
531 frecpoints.SetTime1stC((besttimeC * channelWidth)- 1000.*fLatencyHPTDC +1000.*fLatencyL1C - 1000.*fGRPdelays - fTimeMeanShift[2]);
532 else
533 {
534 frecpoints.SetTimeBestC((besttimeC_best * channelWidth ));
535 frecpoints.SetTime1stC((besttimeC * channelWidth - fTimeMeanShift[2]));
536 }
537 }
538 AliDebug(5,Form("1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f ",
539 besttimeA, besttimeA_best,
540 besttimeC, besttimeC_best) );
541 AliDebug(5,Form("fRecPoints::: 1stimeA %f , besttimeA %f 1sttimeC %f besttimeC %f shiftA %f shiftC %f ",
542 frecpoints.Get1stTimeA(), frecpoints.GetBestTimeA(),
543 frecpoints.Get1stTimeC(), frecpoints.GetBestTimeC(),
544 fTimeMeanShift[1],fTimeMeanShift[2] ) );
545 if( besttimeC < 999999 && besttimeA < 999999) {
c8d939c8 546 if(equalize ==0 )
597a8434 547 timeclock = (channelWidth*(besttimeC + besttimeA)/2.- 1000.*fLatencyHPTDC +1000.*fLatencyL1 - 1000.*fGRPdelays - fTimeMeanShift[0]);
c8d939c8 548 else
549 {
550 timeclock = channelWidth * Float_t( besttimeA+besttimeC)/2. - fTimeMeanShift[0];
551 meanTime = channelWidth * Float_t(besttimeA_best + besttimeC_best )/2.;
552 }
612737bb 553 timeDiff = ( besttimeA - besttimeC)* 0.001* channelWidth ;
9480f05f 554 vertex = meanVertex - c*(timeDiff)/2. ; //+ (fdZonA - fdZonC)/2;
ce50812a 555 }
556
776de217 557 } //if phys event
ce50812a 558 AliDebug(1,Form(" timeDiff %f #channel, meanTime %f #ps, TOFmean%f vertex %f cm meanVertex %f \n",timeDiff, meanTime,timeclock, vertex,meanVertex));
559 frecpoints.SetT0clock(timeclock);
560 frecpoints.SetVertex(vertex);
561 frecpoints.SetMeanTime(meanTime);
562 frecpoints.SetOnlineMean(Int_t(onlineMean));
adf36b9d 563
ce50812a 564 // Set triggers
adf36b9d 565 Bool_t tr[5];
92f6bfd3 566 Int_t trchan[5] = {50,51,52,55,56};
567 Float_t lowtr[5] = {meanTVDC-700, meanOrA-700, meanOrC-700, meanOrC-1000, meanOrC-1000 };
568 Float_t hightr[5] = {meanTVDC+700, meanOrA+700, meanOrC+700, meanOrC+1000, meanOrC+1000};
ce50812a 569
82540e64 570 for (Int_t i=0; i<5; i++) tr[i] = false;
adf36b9d 571 for (Int_t itr=0; itr<5; itr++) {
b0e13b29 572 for (Int_t iHit=0; iHit<1; iHit++)
38cbfa7c 573 {
574 Int_t trr=trchan[itr];
92f6bfd3 575 if( allData[trr][iHit] > lowtr[itr] && allData[trr][iHit] < hightr[itr]) tr[itr]=true;
ce50812a 576
577 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 578 }
adf36b9d 579 }
7e08771f 580 frecpoints.SetT0Trig(tr);
b0e13b29 581
ce50812a 582 // all times with amplitude correction
583 Float_t timecent;
b0e13b29 584 for (Int_t iHit=0; iHit<5; iHit++)
585 {
ce50812a 586 timefull = timecent = -9999;
82540e64 587 tvdc = ora = orc = -9999;
b0e13b29 588 if(allData[50][iHit]>0)
589 tvdc = (Float_t(allData[50][iHit]) - meanTVDC) * channelWidth* 0.001;
590 if(allData[51][iHit]>0)
591 ora = (Float_t(allData[51][iHit]) - meanOrA) * channelWidth* 0.001;
592
593 if(allData[52][iHit]>0)
594 orc = (Float_t(allData[52][iHit]) - meanOrC) * channelWidth* 0.001;
595
7e08771f 596 frecpoints.SetOrC( iHit, orc);
597 frecpoints.SetOrA( iHit, ora);
598 frecpoints.SetTVDC( iHit, tvdc);
b0e13b29 599 for (Int_t i0=0; i0<12; i0++) {
ce50812a 600 if (equalize ==0 )
601 timecent = fTime0vertex[i0] + timeDelayCFD[i0];
602 else
603 timecent = fTime0vertex[i0];
b0e13b29 604 timefull = -9999;
605 if(allData[i0+1][iHit]>1)
ce50812a 606 timefull = (Float_t(allData[i0+1][iHit]) - timecent)* channelWidth* 0.001;
7e08771f 607 frecpoints.SetTimeFull(i0, iHit,timefull) ;
82540e64 608 // 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 609
610 }
611
612 for (Int_t i0=12; i0<24; i0++) {
613 timefull = -9999;
ce50812a 614 if (equalize ==0 )
615 timecent = fTime0vertex[i0] + timeDelayCFD[i0];
616 else
617 timecent = fTime0vertex[i0];
94b4f7e2 618 if(allData[i0+45][iHit]>1) {
ce50812a 619 timefull = (Float_t(allData[i0+45][iHit]) - timecent)* channelWidth* 0.001;
94b4f7e2 620 }
82540e64 621 // 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 622 frecpoints.SetTimeFull(i0, iHit, timefull) ;
b0e13b29 623 }
624 }
625
626
669dc07f 627 //Set MPD
993257ce 628 for (Int_t iHit=0; iHit<5; iHit++)
629 {
630 if (allData[54][iHit] > fTime0vertex[0]+2000 &&
631 allData[54][iHit] <fTime0vertex[0]+3000 &&
632 allData[53][iHit] > fTime0vertex[0]+3000) {
633 frecpoints.SetMultA(allData[53][iHit]-allData[54][iHit]);
634 break;
635 }
636 }
637 for (Int_t iHit=0; iHit<5; iHit++)
638 if (allData[106][iHit] > fTime0vertex[0]+2000 &&
639 allData[106][iHit] <fTime0vertex[0]+3000 &&
640 allData[105][iHit] > fTime0vertex[0]+3000)
641 {
642 frecpoints.SetMultC(allData[105][iHit]-allData[106][iHit]);
643 break;
644 }
b0e13b29 645
669dc07f 646 } // if (else )raw data
58bd3a16 647 recTree->Fill();
58bd3a16 648}
adf36b9d 649
650
ce50812a 651//____________________________________________________________
adf36b9d 652
653 void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
654 {
dc7ca31d 655
656 /***************************************************
657 Resonstruct digits to vertex position
658 ****************************************************/
659
dc7ca31d 660 AliDebug(1,Form("Start FillESD T0"));
b0ab3f59 661 if(!pESD) {
662 AliError("No ESD Event");
663 return;
664 }
4cbe597e 665 pESD ->SetT0spread(fTimeSigmaShift);
666
36cde487 667
58bd3a16 668 Float_t channelWidth = fParam->GetChannelWidth() ;
f7c2c2fc 669 Float_t c = 0.0299792458; // cm/ps
adf36b9d 670 Float_t currentVertex=0, shift=0;
291f31a1 671 Int_t ncont=-1;
adf36b9d 672 const AliESDVertex* vertex = pESD->GetPrimaryVertex();
673 if (!vertex) vertex = pESD->GetPrimaryVertexSPD();
674 if (!vertex) vertex = pESD->GetPrimaryVertexTPC();
675 if (!vertex) vertex = pESD->GetVertex();
676
677 if (vertex) {
678 AliDebug(2, Form("Got %s (%s) from ESD: %f",
679 vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
680 currentVertex = vertex->GetZ();
681
682 ncont = vertex->GetNContributors();
291f31a1 683 if(ncont>0 ) {
adf36b9d 684 shift = currentVertex/c;
adf36b9d 685 }
686 }
d76c31f4 687 TTree *treeR = clustersTree;
dc7ca31d 688
7e08771f 689 AliT0RecPoint frecpoints;
690 AliT0RecPoint * pfrecpoints = &frecpoints;
dc7ca31d 691
692 AliDebug(1,Form("Start FillESD T0"));
693 TBranch *brRec = treeR->GetBranch("T0");
694 if (brRec) {
7e08771f 695 brRec->SetAddress(&pfrecpoints);
dc7ca31d 696 }else{
f16935f7 697 AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d 698 return;
699 }
73df58ab 700
701 brRec->GetEntry(0);
702 Double32_t amp[24], time[24], ampQTC[24], timecorr[24];
703 Double32_t* tcorr;
704 for(Int_t i=0; i<24; i++)
705 amp[i]=time[i]=ampQTC[i]=timecorr[i]=0;
706
ce50812a 707 //1st time
73df58ab 708 Double32_t timeClock[3];
7e08771f 709 Double32_t zPosition = frecpoints.GetVertex();
830172c8 710
ce50812a 711 timeClock[0] = frecpoints.GetT0clock() ;
712 timeClock[1] = frecpoints.Get1stTimeA() + shift;
713 timeClock[2] = frecpoints.Get1stTimeC() - shift;
714 //best time
715 Double32_t timemean[3];
716 timemean[0] = frecpoints.GetMeanTime();
717 timemean[1] = frecpoints.GetBestTimeA() + shift;
718 timemean[2] = frecpoints.GetBestTimeC() - shift;
719
720 for(Int_t i=0; i<3; i++) {
721 fESDTZERO->SetT0TOF(i,timeClock[i]); // interaction time (ns)
722 fESDTZERO->SetT0TOFbest(i,timemean[i]); // interaction time (ns)
723 }
73df58ab 724 for ( Int_t i=0; i<24; i++) {
7e08771f 725 time[i] = frecpoints.GetTime(i); // ps to ns
ce50812a 726 if ( time[i] != 0 && time[i]>-9999) {
7e08771f 727 ampQTC[i] = frecpoints.GetAmp(i);
728 amp[i] = frecpoints.AmpLED(i);
612737bb 729 AliDebug(1,Form("T0: %i time %f ampQTC %f ampLED %f \n", i, time[i], ampQTC[i], amp[i]));
85f61e3b 730 }
73df58ab 731 }
ce50812a 732 fESDTZERO->SetT0time(time); // best TOF on each PMT
733 fESDTZERO->SetT0amplitude(ampQTC); // number of particles(MIPs) on each PMT
7e08771f 734 Int_t trig= frecpoints.GetT0Trig();
735 frecpoints.PrintTriggerSignals( trig);
92f6bfd3 736 // printf(" !!!!! FillESD trigger %i \n",trig);
b0e13b29 737 fESDTZERO->SetT0Trig(trig);
b0e13b29 738 fESDTZERO->SetT0zVertex(zPosition); //vertex Z position
669dc07f 739
7e08771f 740 Double32_t multA=frecpoints.GetMultA();
741 Double32_t multC=frecpoints.GetMultC();
b0e13b29 742 fESDTZERO->SetMultA(multA); // for backward compatubility
743 fESDTZERO->SetMultC(multC); // for backward compatubility
744
745
746 for (Int_t iHit =0; iHit<5; iHit++ ) {
82540e64 747 AliDebug(10,Form("FillESD ::: iHit %i tvdc %f orA %f orC %f\n", iHit,
7e08771f 748 frecpoints.GetTVDC(iHit),
749 frecpoints.GetOrA(iHit),
750 frecpoints.GetOrC(iHit) ));
751 fESDTZERO->SetTVDC(iHit,frecpoints.GetTVDC(iHit));
752 fESDTZERO->SetOrA(iHit,frecpoints.GetOrA(iHit));
753 fESDTZERO->SetOrC(iHit,frecpoints.GetOrC(iHit));
b0e13b29 754
755 for (Int_t i0=0; i0<24; i0++) {
7e08771f 756 // if(frecpoints.GetTimeFull(i0,iHit)>0){
757 // printf("FillESD ::: iHit %i cfd %i time %f \n", iHit, i0, frecpoints.GetTimeFull(i0,iHit));
758 fESDTZERO->SetTimeFull(i0, iHit,frecpoints.GetTimeFull(i0,iHit));
94b4f7e2 759 // }
b0e13b29 760
761 }
762 }
73df58ab 763
ce50812a 764 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 765
ea1a8005 766 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
767 // background flags
768 Bool_t background = BackgroundFlag();
769 fESDTZERO->SetBackgroundFlag(background);
770 Bool_t pileup = PileupFlag();
771 fESDTZERO->SetPileupFlag(pileup);
528890e5 772 for (Int_t i=0; i<5; i++) {
7e08771f 773 fESDTZERO->SetPileupTime(i, frecpoints.GetTVDC(i) ) ;
774 // printf("!!!!!! FillESD :: pileup %i %f %f \n", i,fESDTZERO->GetPileupTime(i), frecpoints.GetTVDC(i));
528890e5 775 }
ea1a8005 776 Bool_t sat = SatelliteFlag();
777 fESDTZERO->SetSatelliteFlag(sat);
73df58ab 778
ea1a8005 779 //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
73df58ab 780 if (pESD) {
adf36b9d 781
73df58ab 782 AliESDfriend *fr = (AliESDfriend*)pESD->FindListObject("AliESDfriend");
783 if (fr) {
82540e64 784 AliDebug(10, Form("Writing TZERO friend data to ESD tree"));
73df58ab 785
85f61e3b 786 // if (ncont>2) {
73df58ab 787 tcorr = fESDTZEROfriend->GetT0timeCorr();
788 for ( Int_t i=0; i<24; i++) {
85f61e3b 789 if(i<12 && time[i]>1) timecorr[i] = tcorr[i] - shift/channelWidth;
790 if(i>11 && time[i]>1) timecorr[i] = tcorr[i] + shift/channelWidth;
612737bb 791 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 792 }
73df58ab 793 fESDTZEROfriend->SetT0timeCorr( timecorr) ;
794 fESDTZEROfriend->SetT0ampLEDminCFD(amp);
795 fESDTZEROfriend->SetT0ampQTC(ampQTC);
85f61e3b 796 fr->SetTZEROfriend(fESDTZEROfriend);
797 // }//
58bd3a16 798 }
b0e13b29 799
800 pESD->SetTZEROData(fESDTZERO);
58bd3a16 801 }
73df58ab 802
dc7ca31d 803} // vertex in 3 sigma
804
ea1a8005 805//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
806 //____________________________________________________________
807
808Bool_t AliT0Reconstructor::PileupFlag() const
809{
810 //
811 Bool_t pileup = false;
528890e5 812 Float_t tvdc[5];
ea1a8005 813 for (Int_t ih=0; ih<5; ih++)
814 {
815 tvdc[ih] = fESDTZERO->GetTVDC(ih);
816
528890e5 817 if( tvdc[0] !=0 && tvdc[0]> -10 && tvdc[0]< 10 )
ea1a8005 818 if(ih>0 && tvdc[ih]>20 ) pileup = true;
528890e5 819 if( tvdc[0] >20 || (tvdc[0] < -20 && tvdc[0] > -9000) ) pileup =true;
820 // if (pileup) printf(" !!!!! pile up %i tvdc %f \n",ih, tvdc[ih]);
ea1a8005 821 }
822
823
824 return pileup;
825
826}
827
828 //____________________________________________________________
829
830Bool_t AliT0Reconstructor::BackgroundFlag() const
831{
9120829d 832
ea1a8005 833 Bool_t background = false;
9120829d 834 /*
835 Float_t orA = fESDTZERO->GetOrA(0);
836 Float_t orC = fESDTZERO->GetOrC(0);
837 Float_t tvdc = fESDTZERO->GetTVDC(ih);
838
839 if ( (orA > -5 && orA <5) && (orC > -5 && orC <5) && (tvdc < -5 || tvdc > 5)) {
840 background = true;
841 // printf(" orA %f orC %f tvdc %f\n", orA, orC, tvdc);
842 }
843 */
ea1a8005 844 return background;
845
846
847}
848
849
850 //____________________________________________________________
851
852Bool_t AliT0Reconstructor::SatelliteFlag() const
853{
854
ce50812a 855 Float_t satelliteLow = GetRecoParam() -> GetLowSatelliteThreshold();
856 Float_t satelliteHigh = GetRecoParam() -> GetHighSatelliteThreshold();
ea1a8005 857 Bool_t satellite = false;
858 for (Int_t i0=0; i0<24; i0++) {
859 Float_t timefull = fESDTZERO -> GetTimeFull(i0,0);
ce50812a 860 if( timefull > satelliteLow && timefull < satelliteHigh) satellite=true;
ea1a8005 861 }
862
863 return satellite;
864
865}