1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 /*********************************************************************
18 * T0 reconstruction and filling ESD
19 * - reconstruct mean time (interation time)
22 ********************************************************************/
24 #include <AliESDEvent.h>
26 #include "AliT0RecPoint.h"
27 #include "AliRawReader.h"
28 #include "AliT0RawReader.h"
29 #include "AliT0digit.h"
30 #include "AliT0Reconstructor.h"
31 #include "AliT0Parameters.h"
32 #include "AliT0Calibrator.h"
37 //#include <iostream.h>
39 ClassImp(AliT0Reconstructor)
41 AliT0Reconstructor:: AliT0Reconstructor(): AliReconstructor(),
50 AliDebug(1,"Start reconstructor ");
52 fParam = AliT0Parameters::Instance();
54 for (Int_t i=0; i<24; i++){
55 TGraph* gr = fParam ->GetAmpLEDRec(i);
56 fAmpLEDrec.AddAtAndExpand(gr,i) ;
58 fdZonC = TMath::Abs(fParam->GetZPositionShift("T0/C/PMT1"));
59 fdZonA = TMath::Abs(fParam->GetZPositionShift("T0/A/PMT15"));
62 //____________________________________________________________________
64 AliT0Reconstructor::AliT0Reconstructor(const AliT0Reconstructor &r):
73 // AliT0Reconstructor copy constructor
76 ((AliT0Reconstructor &) r).Copy(*this);
80 //_____________________________________________________________________________
81 AliT0Reconstructor &AliT0Reconstructor::operator=(const AliT0Reconstructor &r)
84 // Assignment operator
87 if (this != &r) ((AliT0Reconstructor &) r).Copy(*this);
92 //_____________________________________________________________________________
94 void AliT0Reconstructor::Reconstruct(TTree*digitsTree, TTree*clustersTree) const
97 // T0 digits reconstruction
101 TArrayI * timeCFD = new TArrayI(24);
102 TArrayI * timeLED = new TArrayI(24);
103 TArrayI * chargeQT0 = new TArrayI(24);
104 TArrayI * chargeQT1 = new TArrayI(24);
106 AliT0Calibrator *calib=new AliT0Calibrator();
108 // Int_t mV2Mip = param->GetmV2Mip();
109 //mV2Mip = param->GetmV2Mip();
110 Float_t channelWidth = fParam->GetChannelWidth() ;
111 Int_t meanT0 = fParam->GetMeanT0();
113 AliDebug(1,Form("Start DIGITS reconstruction "));
115 TBranch *brDigits=digitsTree->GetBranch("T0");
116 AliT0digit *fDigits = new AliT0digit() ;
118 brDigits->SetAddress(&fDigits);
120 AliError(Form("EXEC Branch T0 digits not found"));
124 digitsTree->GetEvent(0);
125 digitsTree->GetEntry(0);
126 brDigits->GetEntry(0);
127 fDigits->GetTimeCFD(*timeCFD);
128 fDigits->GetTimeLED(*timeLED);
129 fDigits->GetQT0(*chargeQT0);
130 fDigits->GetQT1(*chargeQT1);
133 Float_t besttimeA=999999;
134 Float_t besttimeC=999999;
135 Int_t pmtBestA=99999;
136 Int_t pmtBestC=99999;
137 Float_t timeDiff=999999, meanTime=0;
141 AliT0RecPoint* frecpoints= new AliT0RecPoint ();
142 clustersTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
144 Float_t time[24], adc[24];
145 for (Int_t ipmt=0; ipmt<24; ipmt++) {
146 if(timeCFD->At(ipmt)>0 ){
147 Double_t qt0 = Double_t(chargeQT0->At(ipmt));
148 Double_t qt1 = Double_t(chargeQT1->At(ipmt));
149 if((qt1-qt0)>0) adc[ipmt] = TMath::Exp( Double_t (channelWidth*(qt1-qt0)/1000));
150 time[ipmt] = calib-> WalkCorrection( ipmt, Int_t(qt1) , timeCFD->At(ipmt), "pdc" ) ;
153 Double_t sl = (timeLED->At(ipmt) - time[ipmt])*channelWidth;
154 Double_t qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl/1000.);
155 frecpoints->SetTime(ipmt,time[ipmt]);
156 frecpoints->SetAmp(ipmt,adc[ipmt]);
157 frecpoints->SetAmpLED(ipmt,qt);
165 for (Int_t ipmt=0; ipmt<12; ipmt++){
166 if(time[ipmt] > 1 ) {
167 if(time[ipmt]<besttimeC){
168 besttimeC=time[ipmt]; //timeC
173 for ( Int_t ipmt=12; ipmt<24; ipmt++){
175 if(time[ipmt]<besttimeA) {
176 besttimeA=time[ipmt]; //timeA
180 if(besttimeA !=999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
181 if( besttimeC != 999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
182 AliDebug(1,Form(" besttimeA %f ps, besttimeC %f ps",besttimeA, besttimeC));
183 Float_t c = 0.0299792; // cm/ps
185 if(besttimeA !=999999 && besttimeC != 999999 ){
186 timeDiff =(besttimeC - besttimeA)*channelWidth;
187 meanTime = (meanT0 - (besttimeA + besttimeC)/2) * channelWidth;
188 vertex = c*(timeDiff)/2. + (fdZonA - fdZonC)/2; //-(lenr-lenl))/2;
189 AliDebug(1,Form(" timeDiff %f ps, meanTime %f ps, vertex %f cm",timeDiff, meanTime,vertex ));
190 frecpoints->SetVertex(vertex);
191 frecpoints->SetMeanTime(Int_t(meanTime));
194 //time in each channel as time[ipmt]-MeanTimeinThisChannel(with vertex=0)
195 for (Int_t ipmt=0; ipmt<24; ipmt++) {
197 // time[ipmt] = (time[ipmt] - fTime0vertex[ipmt])*channelWidth;
198 time[ipmt] = time[ipmt] * channelWidth;
199 frecpoints->SetTime(ipmt,time[ipmt]);
202 clustersTree->Fill();
211 //_______________________________________________________________________
213 void AliT0Reconstructor::Reconstruct(AliRawReader* rawReader, TTree*recTree) const
216 // T0RecPoint writing
218 //Q->T-> coefficients !!!! should be measured!!!
219 Int_t allData[110][5];
221 Int_t timeCFD[24], timeLED[24], chargeQT0[24], chargeQT1[24];
222 TString option = GetOption();
223 AliDebug(1,Form("Option: %s\n", option.Data()));
225 for (Int_t i0=0; i0<105; i0++)
227 for (Int_t j0=0; j0<5; j0++) allData[i0][j0]=0;
231 AliDebug(10," before read data ");
232 AliT0RawReader myrawreader(rawReader);
233 if (!myrawreader.Next())
234 AliDebug(1,"No T0 raw data found!");
236 for (Int_t i=0; i<105; i++) {
237 for (Int_t iHit=0; iHit<5; iHit++)
239 allData[i][iHit] = myrawreader.GetData(i,iHit);
242 AliT0Calibrator *calib = new AliT0Calibrator();
244 // Int_t mV2Mip = param->GetmV2Mip();
245 //mV2Mip = param->GetmV2Mip();
246 Float_t channelWidth = fParam->GetChannelWidth() ;
248 Int_t meanT0 = fParam->GetMeanT0();
250 for (Int_t in=0; in<24; in++)
253 timeLED[in] = allData[in+1][0];
254 timeCFD[in] = allData[in+25][0];
255 chargeQT1[in] = allData[in+57][0];
256 chargeQT0[in] = allData[in+80][0];
260 if(option == "cosmic") {
261 for (Int_t in=0; in<12; in++)
263 timeCFD[in] = allData[in+1][0];
264 timeCFD[in+12] = allData[in+1][0];
265 timeLED[in] = allData[in+12+1][0];
266 timeLED[in+12] = allData[in+68+1][0];
269 for (Int_t in=0; in<24; in=in+2)
272 chargeQT1[cc]=allData[in+25][0];
273 chargeQT0[cc]=allData[in+26][0];
275 for (Int_t in=24; in<48; in=in+2)
278 chargeQT1[cc]=allData[in+57][0];
279 chargeQT0[cc]=allData[in+58][0];
283 for (Int_t in=0; in<24; in++)
284 AliDebug(10, Form(" readed Raw %i %i %i %i %i", in, timeLED[in],timeCFD[in],chargeQT0[in],chargeQT1[in]));
287 Float_t besttimeA=9999999;
288 Float_t besttimeC=9999999;
289 Int_t pmtBestA=99999;
290 Int_t pmtBestC=99999;
291 Float_t timeDiff=9999999, meanTime=0;
294 AliT0RecPoint* frecpoints= new AliT0RecPoint ();
296 recTree->Branch( "T0", "AliT0RecPoint" ,&frecpoints, 405,1);
299 Float_t time[24], adc[24];
300 for (Int_t ipmt=0; ipmt<24; ipmt++) {
301 if(timeCFD[ipmt]>0 ){
304 Double_t qt0 = Double_t(chargeQT0[ipmt]);
305 Double_t qt1 = Double_t(chargeQT1[ipmt]);
306 if((qt1-qt0)>0) adc[ipmt] = TMath::Exp( Double_t (channelWidth*(qt1-qt0)/1000));
308 time[ipmt] = calib-> WalkCorrection( ipmt,Int_t(qt1) , timeCFD[ipmt] ) ;
309 Double_t sl = (timeLED[ipmt] - time[ipmt])*channelWidth;
310 if(fAmpLEDrec.At(ipmt))
311 qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl/1000.);
312 frecpoints->SetTime(ipmt,time[ipmt]);
313 frecpoints->SetAmp(ipmt,adc[ipmt]);
314 frecpoints->SetAmpLED(ipmt,qt);
315 AliDebug(1,Form(" QTC %f mv, time in chann %f ",adc[ipmt] ,time[ipmt]));
317 if(option == "cosmic") {
318 Float_t qt0 = Float_t(chargeQT0[ipmt]);
319 Float_t qt1 = Float_t(chargeQT1[ipmt]);
320 if((qt1-qt0)>0) adc[ipmt] = qt1-qt0;
322 time[ipmt] = calib-> WalkCorrection( ipmt, Int_t(adc[ipmt]), timeCFD[ipmt] ) ;
323 Double_t sl = timeLED[ipmt] - time[ipmt];
324 if(fAmpLEDrec.At(ipmt))
325 qt=((TGraph*)fAmpLEDrec.At(ipmt))->Eval(sl/1000.);
326 frecpoints->SetTime(ipmt,time[ipmt]);
327 frecpoints->SetAmp(ipmt,adc[ipmt]);
328 frecpoints->SetAmpLED(ipmt,qt);
329 AliDebug(10,Form(" ipmt %i QTC %i , time in chann %i led %i ",
330 ipmt, Int_t(adc[ipmt]) ,Int_t(time[ipmt]),Int_t( qt)));
340 for (Int_t ipmt=0; ipmt<12; ipmt++){
341 if(time[ipmt] > 1 ) {
342 if(time[ipmt]<besttimeC){
343 besttimeC=time[ipmt]; //timeC
348 for ( Int_t ipmt=12; ipmt<24; ipmt++){
350 if(time[ipmt]<besttimeA) {
351 besttimeA=time[ipmt]; //timeA
355 if(besttimeA !=9999999) frecpoints->SetTimeBestA(Int_t(besttimeA));
356 if( besttimeC != 9999999 ) frecpoints->SetTimeBestC(Int_t(besttimeC));
357 AliDebug(1,Form(" besttimeA %f ps, besttimeC %f ps",besttimeA, besttimeC));
358 Float_t c = 0.0299792; // cm/ps
359 Float_t vertex = 99999;
360 if(besttimeA <9999999 && besttimeC < 9999999 ){
361 timeDiff =( besttimeC - besttimeA) *channelWidth;
363 meanTime = (meanT0 - (besttimeA + besttimeC)/2) * channelWidth;
364 if(option == "cosmic") meanTime = (besttimeA + besttimeC)/2;
365 vertex = c*(timeDiff)/2. + (fdZonA - fdZonC)/2;
366 AliDebug(1,Form(" timeDiff %f ps, meanTime %f ps, vertex %f cm",timeDiff, meanTime,vertex ));
367 frecpoints->SetVertex(vertex);
368 frecpoints->SetMeanTime(Int_t(meanTime));
371 //time in each channel as time[ipmt]-MeanTimeinThisChannel(with vertex=0)
373 for (Int_t ipmt=0; ipmt<24; ipmt++) {
375 time[ipmt] = (time[ipmt] - fTime0vertex[ipmt])*channelWidth;
376 frecpoints->SetTime(ipmt,time[ipmt]);
382 //____________________________________________________________
384 void AliT0Reconstructor::FillESD(TTree */*digitsTree*/, TTree *clustersTree, AliESDEvent *pESD) const
387 /***************************************************
388 Resonstruct digits to vertex position
389 ****************************************************/
391 AliDebug(1,Form("Start FillESD T0"));
393 TTree *treeR = clustersTree;
395 AliT0RecPoint* frecpoints= new AliT0RecPoint ();
397 AliError("Reconstruct Fill ESD >> no recpoints found");
401 AliDebug(1,Form("Start FillESD T0"));
402 TBranch *brRec = treeR->GetBranch("T0");
404 brRec->SetAddress(&frecpoints);
406 AliError(Form("EXEC Branch T0 rec not found"));
411 Float_t amp[24], time[24];
412 Float_t zPosition = frecpoints -> GetVertex();
413 Float_t timeStart = frecpoints -> GetMeanTime() ;
414 for ( Int_t i=0; i<24; i++) {
415 time[i] = Float_t (frecpoints -> GetTime(i)); // ps to ns
416 amp[i] = frecpoints -> GetAmp(i);
418 pESD->SetT0zVertex(zPosition); //vertex Z position
419 pESD->SetT0(timeStart); // interaction time
420 pESD->SetT0time(time); // best TOF on each PMT
421 pESD->SetT0amplitude(amp); // number of particles(MIPs) on each PMT
423 AliDebug(1,Form(" Z position %f cm, T0 %f ps",zPosition , timeStart));
425 } // vertex in 3 sigma