clustering parameter is set to decide the clustering algo
[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
c555f418 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));
135 // time[ipmt] = fCalib-> WalkCorrection( ipmt, Int_t(sl), timeCFD[ipmt],"cosmic" ) ;
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 }
246
c555f418 247 Int_t ref = allData[refPoint][0]-5000.;
bce12dc5 248 Float_t channelWidth = fParam->GetChannelWidth() ;
249
250 // Int_t meanT0 = fParam->GetMeanT0();
d0bcd1fb 251
bce12dc5 252
253 for (Int_t in=0; in<12; in++)
254 {
255 timeCFD[in] = allData[in+1][0] ;
256 timeCFD[in+12] = allData[in+56+1][0] ;
257 timeLED[in] = allData[in+12+1][0] ;
258 timeLED[in+12] = allData[in+68+1][0] ;
d0bcd1fb 259 AliDebug(10, Form(" readed i %i cfdC %i cfdA %i ledC %i ledA%i ",
260 in, timeCFD[in],timeCFD[in+12],timeLED[in],
261 timeLED[in+12]));
bce12dc5 262 }
263
264 for (Int_t in=0; in<12; in++)
265 {
d0bcd1fb 266 chargeQT0[in]=allData[2*in+25][0];
267 chargeQT1[in]=allData[2*in+26][0];
bce12dc5 268 }
269
270 for (Int_t in=12; in<24; in++)
271 {
d0bcd1fb 272 chargeQT0[in]=allData[2*in+57][0];
273 chargeQT1[in]=allData[2*in+58][0];
bce12dc5 274 }
275
d0bcd1fb 276 // } //cosmic with physics event
2e6a5ee0 277 for (Int_t in=0; in<24; in++)
c883fdf2 278 AliDebug(10, Form(" readed Raw %i %i %i %i %i",
279 in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
776de217 280 onlineMean = allData[49][0];
2e6a5ee0 281
12e9daf9 282 Float_t time[24], adc[24];
2e6a5ee0 283 for (Int_t ipmt=0; ipmt<24; ipmt++) {
284 if(timeCFD[ipmt]>0 && timeLED[ipmt]>0){
bce12dc5 285 //for simulated data
d0bcd1fb 286 //for physics data
287 if(( chargeQT1[ipmt] - chargeQT0[ipmt])>0)
288 adc[ipmt] = chargeQT1[ipmt] - chargeQT0[ipmt];
289 else
290 adc[ipmt] = 0;
bce12dc5 291
539b9cb9 292
40043664 293 // time[ipmt] = fCalib-> WalkCorrection(refAmp, ipmt, adc[ipmt], timeCFD[ipmt] ) ;
bce12dc5 294
d0bcd1fb 295 Double_t sl = timeLED[ipmt] - timeCFD[ipmt];
40043664 296 time[ipmt] = fCalib-> WalkCorrection( refAmp,ipmt, Int_t(sl), timeCFD[ipmt] ) ;
d0bcd1fb 297 AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i (led-cfd) %i ",
298 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( sl)));
299 Double_t ampMip =( (TGraph*)fAmpLED.At(ipmt))->Eval(sl);
300 Double_t qtMip = ((TGraph*)fQTC.At(ipmt))->Eval(adc[ipmt]);
301 AliDebug(10,Form(" Amlitude in MIPS LED %f ; QTC %f; in channels %i\n ",ampMip,qtMip, adc[ipmt]));
bce12dc5 302
d0bcd1fb 303 frecpoints->SetTime(ipmt, Float_t(time[ipmt]) );
304 frecpoints->SetAmp(ipmt, Float_t( ampMip)); //for cosmic &pp beam
305 frecpoints->SetAmpLED(ipmt, Float_t(qtMip));
bce12dc5 306
2e6a5ee0 307 }
308 else {
309 time[ipmt] = 0;
310 adc[ipmt] = 0;
311 }
312 }
313
314 for (Int_t ipmt=0; ipmt<12; ipmt++){
315 if(time[ipmt] > 1 ) {
316 if(time[ipmt]<besttimeC){
317 besttimeC=time[ipmt]; //timeC
318 pmtBestC=ipmt;
319 }
320 }
321 }
322 for ( Int_t ipmt=12; ipmt<24; ipmt++){
323 if(time[ipmt] > 1) {
324 if(time[ipmt]<besttimeA) {
325 besttimeA=time[ipmt]; //timeA
326 pmtBestA=ipmt;}
327 }
328 }
329 if(besttimeA !=9999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
330 if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
539b9cb9 331 AliDebug(5,Form(" pmtA %i besttimeA %f ps, pmtC %i besttimeC %f #channel",
d0bcd1fb 332 pmtBestA,besttimeA, pmtBestC, besttimeC));
776de217 333 if(besttimeA <9999999 && besttimeC < 9999999 ){
9b83615d 334 timeDiff = ( besttimeA - besttimeC) *channelWidth;
d0bcd1fb 335 meanTime = Float_t((besttimeA + besttimeC)/2.);
776de217 336 timeclock = Float_t(meanTime - ref);
d0bcd1fb 337 vertex = meanVertex - c*(timeDiff)/2.; //+ (fdZonA - fdZonC)/2;
2e6a5ee0 338 }
776de217 339 } //if phys event
340 frecpoints->SetT0clock(timeclock);
341 frecpoints->SetVertex(vertex);
342 frecpoints->SetMeanTime(Int_t(meanTime));
343 frecpoints->SetOnlineMean(Int_t(onlineMean));
344 AliDebug(5,Form(" timeDiff %f #channel, meanTime %f #channel, vertex %f cm meanVertex %f online mean %i ",timeDiff, meanTime,vertex,meanVertex, onlineMean));
345
346
bce12dc5 347 } // if (else )raw data
348 recTree->Fill();
349 if(frecpoints) delete frecpoints;
c41ceaac 350}
bce12dc5 351
352
c41ceaac 353//____________________________________________________________
354
d76c31f4 355void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
dc7ca31d 356{
357
358 /***************************************************
359 Resonstruct digits to vertex position
360 ****************************************************/
361
dc7ca31d 362 AliDebug(1,Form("Start FillESD T0"));
363
d76c31f4 364 TTree *treeR = clustersTree;
dc7ca31d 365
366 AliT0RecPoint* frecpoints= new AliT0RecPoint ();
367 if (!frecpoints) {
368 AliError("Reconstruct Fill ESD >> no recpoints found");
369 return;
370 }
371
372 AliDebug(1,Form("Start FillESD T0"));
373 TBranch *brRec = treeR->GetBranch("T0");
374 if (brRec) {
375 brRec->SetAddress(&frecpoints);
376 }else{
f16935f7 377 AliError(Form("EXEC Branch T0 rec not found"));
dc7ca31d 378 return;
379 }
380
381 brRec->GetEntry(0);
f16935f7 382 Float_t amp[24], time[24];
383 Float_t zPosition = frecpoints -> GetVertex();
384 Float_t timeStart = frecpoints -> GetMeanTime() ;
776de217 385 Float_t timeClock = frecpoints -> GetT0clock() ;
386 for ( Int_t i=0; i<24; i++) {
94c27e4f 387 time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns
dc7ca31d 388 amp[i] = frecpoints -> GetAmp(i);
389 }
f16935f7 390 pESD->SetT0zVertex(zPosition); //vertex Z position
dc7ca31d 391 pESD->SetT0(timeStart); // interaction time
776de217 392 pESD->SetT0clock(timeClock); // interaction time with ref.point(spectrum)
dc7ca31d 393 pESD->SetT0time(time); // best TOF on each PMT
394 pESD->SetT0amplitude(amp); // number of particles(MIPs) on each PMT
5325480c 395
f16935f7 396 AliDebug(1,Form(" Z position %f cm, T0 %f ps",zPosition , timeStart));
dc7ca31d 397
dc7ca31d 398} // vertex in 3 sigma
399
400
401
402
403
404