Bug fixes.
[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"
dc7ca31d 33
34#include <TArrayI.h>
35#include <TGraph.h>
aad72f45 36#include <TMath.h>
b09247a2 37#include <Riostream.h>
dc7ca31d 38
39ClassImp(AliT0Reconstructor)
40
c41ceaac 41 AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
f16935f7 42 fdZonA(0),
43 fdZonC(0),
44 fZposition(0),
45 fParam(NULL),
2e6a5ee0 46 fAmpLEDrec(),
c883fdf2 47 fQTC(0),
48 fAmpLED(0),
2e6a5ee0 49 fCalib()
e0bba6cc 50{
72e48d95 51 //constructor
52
74adb36a 53 fParam = AliT0Parameters::Instance();
54 fParam->Init();
c883fdf2 55
74adb36a 56 for (Int_t i=0; i<24; i++){
2e6a5ee0 57 TGraph* gr = fParam ->GetAmpLEDRec(i);
29ed1d0e 58 if (gr) fAmpLEDrec.AddAtAndExpand(gr,i) ;
c883fdf2 59 TGraph* gr1 = fParam ->GetAmpLED(i);
60 if (gr1) fAmpLED.AddAtAndExpand(gr1,i) ;
61 TGraph* gr2 = fParam ->GetQTC(i);
539b9cb9 62 if (gr2) fQTC.AddAtAndExpand(gr2,i) ;
c883fdf2 63 }
539b9cb9 64
29ed1d0e 65
f16935f7 66 fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
67 fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
539b9cb9 68
12e9daf9 69 fCalib = new AliT0Calibrator();
70
c41ceaac 71}
c41ceaac 72
73//_____________________________________________________________________________
dc7ca31d 74void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
94c27e4f 75
dc7ca31d 76{
94c27e4f 77 // T0 digits reconstruction
c555f418 78 Int_t refAmp = GetRecoParam()->GetRefAmp();
776de217 79
c41ceaac 80 TArrayI * timeCFD = new TArrayI(24);
81 TArrayI * timeLED = new TArrayI(24);
82 TArrayI * chargeQT0 = new TArrayI(24);
83 TArrayI * chargeQT1 = new TArrayI(24);
74adb36a 84
d0bcd1fb 85
8955c6b4 86 Float_t channelWidth = fParam->GetChannelWidth() ;
b95e8d87 87 Float_t meanVertex = fParam->GetMeanVertex();
776de217 88 Float_t c = 0.0299792; // cm/ps
89 Float_t vertex = 9999999;
90 Int_t timeDiff=999999, meanTime=999999, timeclock=999999;
91
94c27e4f 92
dc7ca31d 93 AliDebug(1,Form("Start DIGITS reconstruction "));
94c27e4f 94
2e6a5ee0 95
d0bcd1fb 96 TBranch *brDigits=digitsTree->GetBranch("T0");
e0bba6cc 97 AliT0digit *fDigits = new AliT0digit() ;
dc7ca31d 98 if (brDigits) {
99 brDigits->SetAddress(&fDigits);
100 }else{
f16935f7 101 AliError(Form("EXEC Branch T0 digits not found"));
102 return;
dc7ca31d 103 }
e0bba6cc 104
c41ceaac 105 digitsTree->GetEvent(0);
106 digitsTree->GetEntry(0);
107 brDigits->GetEntry(0);
108 fDigits->GetTimeCFD(*timeCFD);
109 fDigits->GetTimeLED(*timeLED);
110 fDigits->GetQT0(*chargeQT0);
111 fDigits->GetQT1(*chargeQT1);
446d6ec4 112 Int_t onlineMean = fDigits->MeanTime();
776de217 113 Int_t ref = fDigits->RefPoint();
c883fdf2 114
c41ceaac 115
116 Float_t besttimeA=999999;
117 Float_t besttimeC=999999;
118 Int_t pmtBestA=99999;
119 Int_t pmtBestC=99999;
dc7ca31d 120
94c27e4f 121 AliT0RecPoint* frecpoints= new AliT0RecPoint ();
122 clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
123
b95e8d87 124 Float_t time[24], adc[24];
dc7ca31d 125 for (Int_t ipmt=0; ipmt<24; ipmt++) {
c41ceaac 126 if(timeCFD->At(ipmt)>0 ){
d0bcd1fb 127 if(( chargeQT1->At(ipmt) - chargeQT0->At(ipmt))>0)
128 adc[ipmt] = chargeQT1->At(ipmt) - chargeQT0->At(ipmt);
a7027400 129 else
d0bcd1fb 130 adc[ipmt] = 0;
131
541b42c4 132 // time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD->At(ipmt)) ;
d0bcd1fb 133
134 Double_t sl = Double_t(timeLED->At(ipmt) - timeCFD->At(ipmt));
541b42c4 135 time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD->At(ipmt) ) ;
d0bcd1fb 136 AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
137 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
776de217 138
d0bcd1fb 139 Double_t ampMip =((TGraph*)fAmpLED.At(ipmt))->Eval(sl);
140 Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
141 AliDebug(10,Form(" Amlitude in MIPS LED %f , QTC %f in channels %i\n ",ampMip,qtMip, adc[ipmt]));
142
143 frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
144 frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam
145 frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
146
dc7ca31d 147 }
148 else {
149 time[ipmt] = 0;
150 adc[ipmt] = 0;
151 }
152 }
94c27e4f 153
dc7ca31d 154 for (Int_t ipmt=0; ipmt<12; ipmt++){
155 if(time[ipmt] > 1 ) {
c41ceaac 156 if(time[ipmt]<besttimeC){
157 besttimeC=time[ipmt]; //timeC
158 pmtBestC=ipmt;
dc7ca31d 159 }
160 }
161 }
162 for ( Int_t ipmt=12; ipmt<24; ipmt++){
163 if(time[ipmt] > 1) {
c41ceaac 164 if(time[ipmt]<besttimeA) {
165 besttimeA=time[ipmt]; //timeA
166 pmtBestA=ipmt;}
dc7ca31d 167 }
168 }
c41ceaac 169 if(besttimeA !=999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
170 if( besttimeC != 999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
d0bcd1fb 171 AliDebug(10,Form(" besttimeA %f ch, besttimeC %f ch",besttimeA, besttimeC));
776de217 172 if(besttimeA !=999999 && besttimeC != 999999 ){
9b83615d 173 // timeDiff = (besttimeC - besttimeA)*channelWidth;
174 timeDiff = (besttimeA - besttimeC)*channelWidth;
d0bcd1fb 175 meanTime = Float_t((besttimeA + besttimeC)/2);// * channelWidth);
776de217 176 timeclock = Float_t(meanTime - ref);
d0bcd1fb 177 vertex = meanVertex - c*(timeDiff)/2.;// + (fdZonA - fdZonC)/2;
776de217 178 }
dc7ca31d 179 frecpoints->SetVertex(vertex);
776de217 180 frecpoints->SetMeanTime(meanTime);
181 frecpoints->SetT0clock(timeclock);
446d6ec4 182 //online mean
d0bcd1fb 183 frecpoints->SetOnlineMean(Int_t(onlineMean));
776de217 184 AliDebug(10,Form(" timeDiff %i #channel, meanTime %i #channel, vertex %f cm online mean %i timeclock %i ps",timeDiff, meanTime,vertex, Int_t(onlineMean), timeclock));
dc7ca31d 185
776de217 186 // }
b95e8d87 187
dc7ca31d 188 clustersTree->Fill();
bd375212 189
190 delete timeCFD;
191 delete timeLED;
192 delete chargeQT0;
193 delete chargeQT1;
dc7ca31d 194}
195
196
c41ceaac 197//_______________________________________________________________________
198
199void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
200{
94c27e4f 201 // T0 raw ->
539b9cb9 202 //
203 // reference amplitude and time ref. point from reco param
204
c555f418 205 Int_t refAmp = GetRecoParam()->GetRefAmp();
206 Int_t refPoint = GetRecoParam()->GetRefPoint();
539b9cb9 207
e8ed1cd0 208 Int_t allData[110][5];
2e6a5ee0 209
e8ed1cd0 210 Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
776de217 211 Int_t timeDiff=999999, meanTime=999999, timeclock=999999;
212 Float_t c = 0.0299792458; // cm/ps
213 Float_t vertex = 9999999;
214 Int_t onlineMean=0;
bce12dc5 215 for (Int_t i0=0; i0<105; i0++)
216 {
217 for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
218 }
2e6a5ee0 219
bce12dc5 220 Float_t besttimeA=9999999;
221 Float_t besttimeC=9999999;
222 Int_t pmtBestA=99999;
223 Int_t pmtBestC=99999;
b95e8d87 224 Float_t meanVertex = fParam->GetMeanVertex();
29a60970 225
bce12dc5 226 AliT0RecPoint* frecpoints= new AliT0RecPoint ();
227
228 recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
2e6a5ee0 229
bce12dc5 230 AliDebug(10," before read data ");
231 AliT0RawReader myrawreader(rawReader);
776de217 232
233 UInt_t type =rawReader->GetType();
234
bce12dc5 235 if (!myrawreader.Next())
236 AliDebug(1,Form(" no raw data found!!"));
237 else
238 {
776de217 239 if(type == 7) { //only physics
bce12dc5 240 for (Int_t i=0; i<105; i++) {
241 for (Int_t iHit=0; iHit<5; iHit++)
242 {
243 allData[i][iHit] = myrawreader.GetData(i,iHit);
244 }
245 }
541b42c4 246 Int_t ref=0;
247 if (refPoint>0)
248 ref = allData[refPoint][0]-5000.;
249
bce12dc5 250 Float_t channelWidth = fParam->GetChannelWidth() ;
251
252 // Int_t meanT0 = fParam->GetMeanT0();
d0bcd1fb 253
bce12dc5 254
255 for (Int_t in=0; in<12; in++)
256 {
257 timeCFD[in] = allData[in+1][0] ;
258 timeCFD[in+12] = allData[in+56+1][0] ;
259 timeLED[in] = allData[in+12+1][0] ;
260 timeLED[in+12] = allData[in+68+1][0] ;
d0bcd1fb 261 AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
262 in, timeCFD[in],timeCFD[in+12],timeLED[in],
263 timeLED[in+12]));
bce12dc5 264 }
265
266 for (Int_t in=0; in<12; in++)
267 {
d0bcd1fb 268 chargeQT0[in]=allData[2*in+25][0];
269 chargeQT1[in]=allData[2*in+26][0];
bce12dc5 270 }
271
272 for (Int_t in=12; in<24; in++)
273 {
d0bcd1fb 274 chargeQT0[in]=allData[2*in+57][0];
275 chargeQT1[in]=allData[2*in+58][0];
bce12dc5 276 }
277
d0bcd1fb 278 // } //cosmic with physics event
2e6a5ee0 279 for (Int_t in=0; in<24; in++)
c883fdf2 280 AliDebug(10, Form(" readed Raw %i %i %i %i %i",
281 in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
776de217 282 onlineMean = allData[49][0];
2e6a5ee0 283
12e9daf9 284 Float_t time[24], adc[24];
2e6a5ee0 285 for (Int_t ipmt=0; ipmt<24; ipmt++) {
286 if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
bce12dc5 287 //for simulated data
d0bcd1fb 288 //for physics data
289 if(( chargeQT1[ipmt] - chargeQT0[ipmt])>0)
541b42c4 290 adc[ipmt] = chargeQT0[ipmt] - chargeQT1[ipmt];
d0bcd1fb 291 else
292 adc[ipmt] = 0;
bce12dc5 293
539b9cb9 294
40043664 295 // time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ;
bce12dc5 296
d0bcd1fb 297 Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
40043664 298 time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
d0bcd1fb 299 AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
300 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
301 Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
302 Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
303 AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %i\n ",ampMip,qtMip, adc[ipmt]));
bce12dc5 304
d0bcd1fb 305 frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
541b42c4 306 frecpoints->SetAmpLED(ipmt, Float_t( qtMip)); //for cosmic &pp beam
307 frecpoints->SetAmp(ipmt, Float_t(ampMip));
bce12dc5 308
2e6a5ee0 309 }
310 else {
311 time[ipmt] = 0;
312 adc[ipmt] = 0;
313 }
314 }
315
316 for (Int_t ipmt=0; ipmt<12; ipmt++){
317 if(time[ipmt] > 1 ) {
318 if(time[ipmt]<besttimeC){
319 besttimeC=time[ipmt]; //timeC
320 pmtBestC=ipmt;
321 }
322 }
323 }
324 for ( Int_t ipmt=12; ipmt<24; ipmt++){
325 if(time[ipmt] > 1) {
326 if(time[ipmt]<besttimeA) {
327 besttimeA=time[ipmt]; //timeA
328 pmtBestA=ipmt;}
329 }
330 }
331 if(besttimeA !=9999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
332 if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
539b9cb9 333 AliDebug(5,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f #channel",
d0bcd1fb 334 pmtBestA,besttimeA, pmtBestC, besttimeC));
776de217 335 if(besttimeA <9999999 && besttimeC < 9999999 ){
9b83615d 336 timeDiff = ( besttimeA - besttimeC) *channelWidth;
541b42c4 337 meanTime = Float_t(((besttimeA-ref) + (besttimeC-ref))/2.);
338 timeclock = Float_t(meanTime);
d0bcd1fb 339 vertex = meanVertex - c*(timeDiff)/2.; //+ (fdZonA - fdZonC)/2;
2e6a5ee0 340 }
776de217 341 } //if phys event
342 frecpoints->SetT0clock(timeclock);
343 frecpoints->SetVertex(vertex);
344 frecpoints->SetMeanTime(Int_t(meanTime));
345 frecpoints->SetOnlineMean(Int_t(onlineMean));
346 AliDebug(5,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm meanVertex %f online mean %i ",timeDiff, meanTime,vertex,meanVertex, onlineMean));
347
348
bce12dc5 349 } // if (else )raw data
350 recTree->Fill();
351 if(frecpoints) delete frecpoints;
c41ceaac 352}
bce12dc5 353
354
c41ceaac 355//____________________________________________________________
356
d76c31f4 357void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
dc7ca31d 358{
359
360 /***************************************************
361 Resonstruct digits to vertex position
362 ****************************************************/
363
dc7ca31d 364 AliDebug(1,Form("Start FillESD T0"));
365
d76c31f4 366 TTree *treeR = clustersTree;
dc7ca31d 367
368 AliT0RecPoint* frecpoints= new AliT0RecPoint ();
369 if (!frecpoints) {
370 AliError("Reconstruct Fill ESD >> no recpoints found");
371 return;
372 }
373
374 AliDebug(1,Form("Start FillESD T0"));
375 TBranch *brRec = treeR->GetBranch("T0");
376 if (brRec) {
377 brRec->SetAddress(&frecpoints);
378 }else{
f16935f7 379 AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d 380 return;
381 }
382
383 brRec->GetEntry(0);
f16935f7 384 Float_t amp[24], time[24];
385 Float_t zPosition = frecpoints -> GetVertex();
386 Float_t timeStart = frecpoints -> GetMeanTime() ;
776de217 387 Float_t timeClock = frecpoints -> GetT0clock() ;
388 for ( Int_t i=0; i<24; i++) {
94c27e4f 389 time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns
dc7ca31d 390 amp[i] = frecpoints -> GetAmp(i);
391 }
f16935f7 392 pESD->SetT0zVertex(zPosition); //vertex Z position
dc7ca31d 393 pESD->SetT0(timeStart); // interaction time
776de217 394 pESD->SetT0clock(timeClock); // interaction time with ref.point(spectrum)
dc7ca31d 395 pESD->SetT0time(time); // best TOF on each PMT
396 pESD->SetT0amplitude(amp); // number of particles(MIPs) on each PMT
5325480c 397
f16935f7 398 AliDebug(1,Form(" Z position %f cm, T0 %f ps",zPosition , timeStart));
dc7ca31d 399
dc7ca31d 400} // vertex in 3 sigma
401
402
403
404
405
406