Line data Source code
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: AliADDigitizer.cxx $ */
17 :
18 : ///_________________________________________________________________________
19 : ///
20 : /// This class constructs Digits out of Hits
21 : ///
22 : ///
23 :
24 : // --- Standard library ---
25 :
26 : // --- ROOT system ---
27 : #include <TMath.h>
28 : #include <TTree.h>
29 : #include <TMap.h>
30 : #include <TGeoManager.h>
31 : #include <TGeoPhysicalNode.h>
32 : #include <AliGeomManager.h>
33 : #include <TRandom.h>
34 : #include <TF1.h>
35 : #include <TH1F.h>
36 : #include <TH2F.h>
37 : #include <TCanvas.h>
38 :
39 : // --- AliRoot header files ---
40 : #include "AliRun.h"
41 : #include "AliDetector.h"
42 : #include "AliAD.h"
43 : #include "AliADhit.h"
44 : #include "AliADConst.h"
45 : #include "AliRunLoader.h"
46 : #include "AliLoader.h"
47 : #include "AliGRPObject.h"
48 : #include "AliDigitizationInput.h"
49 : #include "AliCDBManager.h"
50 : #include "AliCDBStorage.h"
51 : #include "AliCDBEntry.h"
52 : #include "AliADCalibData.h"
53 : #include "AliCTPTimeParams.h"
54 : #include "AliLHCClockPhase.h"
55 : #include "TSpline.h"
56 : #include "AliADdigit.h"
57 : #include "AliADDigitizer.h"
58 : #include "AliADSDigit.h"
59 : #include "AliADTriggerSimulator.h"
60 : #include "AliLog.h"
61 :
62 12 : ClassImp(AliADDigitizer)
63 :
64 : //____________________________________________________________________________
65 : AliADDigitizer::AliADDigitizer()
66 0 : :AliDigitizer(),
67 0 : fCalibData(GetCalibData()),
68 0 : fNdigits(0),
69 0 : fDigits(0),
70 0 : fTimeSignalShape(NULL),
71 0 : fChargeSignalShape(NULL),
72 0 : fEvenOrOdd(kFALSE),
73 0 : fTask(kHits2Digits),
74 0 : fAD(NULL)
75 0 : {
76 : // default constructor
77 : // Initialize OCDB and containers used in the digitization
78 :
79 0 : Init();
80 0 : }
81 :
82 : //____________________________________________________________________________
83 : AliADDigitizer::AliADDigitizer(AliAD *AD, DigiTask_t task)
84 0 : :AliDigitizer(),
85 0 : fCalibData(GetCalibData()),
86 0 : fNdigits(0),
87 0 : fDigits(0),
88 0 : fTimeSignalShape(NULL),
89 0 : fChargeSignalShape(NULL),
90 0 : fEvenOrOdd(kFALSE),
91 0 : fTask(task),
92 0 : fAD(AD)
93 0 : {
94 : // constructor
95 : // Initialize OCDB and containers used in the digitization
96 :
97 0 : Init();
98 0 : }
99 :
100 : //____________________________________________________________________________
101 : AliADDigitizer::AliADDigitizer(AliDigitizationInput* digInput)
102 0 : :AliDigitizer(digInput),
103 0 : fCalibData(GetCalibData()),
104 0 : fNdigits(0),
105 0 : fDigits(0),
106 0 : fTimeSignalShape(NULL),
107 0 : fChargeSignalShape(NULL),
108 0 : fEvenOrOdd(kFALSE),
109 0 : fTask(kHits2Digits),
110 0 : fAD(NULL)
111 0 : {
112 : // constructor
113 : // Initialize OCDB and containers used in the digitization
114 :
115 0 : Init();
116 0 : }
117 :
118 : //____________________________________________________________________________
119 : AliADDigitizer::~AliADDigitizer()
120 0 : {
121 : // destructor
122 :
123 0 : if (fDigits) {
124 0 : fDigits->Delete();
125 0 : delete fDigits;
126 0 : fDigits=0;
127 0 : }
128 :
129 0 : if (fTimeSignalShape) {
130 0 : delete fTimeSignalShape;
131 0 : fTimeSignalShape = NULL;
132 0 : }
133 0 : if (fChargeSignalShape) {
134 0 : delete fChargeSignalShape;
135 0 : fChargeSignalShape = NULL;
136 0 : }
137 :
138 0 : for(Int_t i = 0 ; i < 16; ++i) {
139 0 : if (fTime[i]) delete [] fTime[i];
140 : }
141 0 : }
142 :
143 : //____________________________________________________________________________
144 : Bool_t AliADDigitizer::Init()
145 : {
146 : // Initialises the digitizer
147 : // Initialize OCDB and containers used in the digitization
148 :
149 : // check if the digitizer was already initialized
150 0 : if (fTimeSignalShape) return kTRUE;
151 0 : fTimeSignalShape = new TF1("ADTimeSignalShape",this,&AliADDigitizer::TimeSignalShape,0,200,6,"AliADDigitizer","TimeSignalShape");
152 0 : fChargeSignalShape = new TF1("ADChargeSignalShape",this,&AliADDigitizer::ChargeSignalShape,0,300,3,"AliADDigitizer","ChargeSignalShape");
153 0 : fThresholdShape = new TF1("ADThresholdShape",this,&AliADDigitizer::ThresholdShape,0,50,1,"AliADDigitizer","ThresholdShape");
154 :
155 0 : fTimeSignalShape->SetParameters(-1.07335e+00,2.16002e+01,-1.26133e-01,
156 : 1.41619e+00,5.50334e-01,3.86111e-01);
157 :
158 : // Now get the CTP L0->L1 delay
159 0 : AliCDBEntry *entry = AliCDBManager::Instance()->Get("GRP/CTP/CTPtiming");
160 0 : if (!entry) AliFatal("CTP timing parameters are not found in OCDB !");
161 0 : AliCTPTimeParams *ctpParams = (AliCTPTimeParams*)entry->GetObject();
162 0 : Float_t l1Delay = (Float_t)ctpParams->GetDelayL1L0()*25.0;
163 :
164 0 : AliCDBEntry *entry1 = AliCDBManager::Instance()->Get("GRP/CTP/TimeAlign");
165 0 : if (!entry1) AliFatal("CTP time-alignment is not found in OCDB !");
166 0 : AliCTPTimeParams *ctpTimeAlign = (AliCTPTimeParams*)entry1->GetObject();
167 0 : l1Delay += ((Float_t)ctpTimeAlign->GetDelayL1L0()*25.0);
168 :
169 0 : AliCDBEntry *entry2 = AliCDBManager::Instance()->Get("AD/Calib/TimeDelays");
170 0 : if (!entry2) AliFatal("AD time delays are not found in OCDB !");
171 0 : TH1F *TimeDelays = (TH1F*)entry2->GetObject();
172 :
173 0 : AliCDBEntry *entry3 = AliCDBManager::Instance()->Get("GRP/Calib/LHCClockPhase");
174 0 : if (!entry3) AliFatal("LHC clock-phase shift is not found in OCDB !");
175 0 : AliLHCClockPhase *phase = (AliLHCClockPhase*)entry3->GetObject();
176 :
177 : //Get Pulse shape parameters
178 0 : AliCDBEntry *entry4 = AliCDBManager::Instance()->Get("AD/Calib/PulseShapes");
179 0 : if (!entry4) AliFatal("AD pulse shapes are not found in OCDB !");
180 0 : TH2F *PulseShapes = (TH2F*)entry4->GetObject();
181 :
182 : //Time slewing splines
183 0 : GetTimeSlewingSplines();
184 : //Try to extrapolate the splines
185 0 : ExtrapolateSplines();
186 :
187 0 : for(Int_t i = 0 ; i < 16; ++i) {
188 :
189 0 : fCssOffset[i] = PulseShapes->GetBinContent(i+1,1);
190 0 : fCssTau[i] = PulseShapes->GetBinContent(i+1,2);
191 0 : fCssSigma[i] = PulseShapes->GetBinContent(i+1,3);
192 :
193 0 : for(Int_t j = 0; j < kADNClocks; ++j) fAdc[i][j] = 0;
194 0 : fLeadingTime[i] = fTimeWidth[i] = 0;
195 :
196 0 : fPmGain[i] = fCalibData->GetGain(i);
197 :
198 0 : fAdcPedestal[i][0] = fCalibData->GetPedestal(i);
199 0 : fAdcSigma[i][0] = fCalibData->GetSigma(i);
200 0 : fAdcPedestal[i][1] = fCalibData->GetPedestal(i+16);
201 0 : fAdcSigma[i][1] = fCalibData->GetSigma(i+16);
202 :
203 0 : Int_t board = AliADCalibData::GetBoardNumber(i);
204 0 : fNBins[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0+
205 0 : (Float_t)kADMaxTDCWidth*fCalibData->GetWidthResolution(board))/
206 0 : fCalibData->GetTimeResolution(board));
207 0 : fNBinsLT[i] = TMath::Nint(((Float_t)(fCalibData->GetMatchWindow(board)+1)*25.0)/
208 0 : fCalibData->GetTimeResolution(board));
209 0 : fBinSize[i] = fCalibData->GetTimeResolution(board);
210 :
211 0 : fHptdcOffset[i] = (((Float_t)fCalibData->GetRollOver(board)-
212 0 : (Float_t)fCalibData->GetTriggerCountOffset(board))*25.0
213 0 : +fCalibData->GetTimeOffset(i)
214 0 : -l1Delay
215 0 : -phase->GetMeanPhase()
216 0 : -TimeDelays->GetBinContent(i+1)
217 0 : -kADOffset);
218 :
219 0 : fClockOffset[i] = (((Float_t)fCalibData->GetRollOver(board)-
220 0 : (Float_t)fCalibData->GetTriggerCountOffset(board))*25.0
221 0 : +fCalibData->GetTimeOffset(i)
222 0 : -l1Delay
223 0 : -kADOffset);
224 :
225 0 : fTime[i] = new Float_t[fNBins[i]];
226 0 : memset(fTime[i],0,fNBins[i]*sizeof(Float_t));
227 :
228 : //std::cout<<"AD: "<<" fNBins = "<<fNBins[i]<<" fNBinsLT = "<<fNBinsLT[i]<<" fHptdcOffset = "<<fHptdcOffset[i]<<" fClockOffset = "<<fClockOffset[i]<<std::endl;
229 : }
230 :
231 : return kTRUE;
232 :
233 0 : }
234 :
235 : //____________________________________________________________________________
236 : void AliADDigitizer::Digitize(Option_t* /*option*/)
237 : {
238 : // Creates digits from hits
239 0 : fNdigits = 0;
240 :
241 0 : if (fAD && !fDigInput) {
242 0 : AliLoader *loader = fAD->GetLoader();
243 0 : if (!loader) {
244 0 : AliError("Can not get AD Loader via AliAD object!");
245 0 : return;
246 : }
247 0 : AliRunLoader* runLoader = AliRunLoader::Instance();
248 0 : for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); ++iEvent) {
249 0 : runLoader->GetEvent(iEvent);
250 0 : if (fTask == kHits2Digits) {
251 0 : DigitizeHits();
252 0 : DigitizeSDigits();
253 0 : WriteDigits(loader);
254 0 : }
255 : else {
256 : DigitizeHits();
257 0 : WriteSDigits(loader);
258 : }
259 : }
260 0 : }
261 0 : else if (fDigInput) {
262 0 : ReadSDigits();
263 0 : DigitizeSDigits();
264 0 : AliRunLoader *currentLoader = AliRunLoader::GetRunLoader(fDigInput->GetOutputFolderName());
265 0 : AliLoader *loader = currentLoader->GetLoader("ADLoader");
266 0 : if (!loader) {
267 0 : AliError("Cannot get AD Loader via RunDigitizer!");
268 0 : return;
269 : }
270 0 : WriteDigits(loader);
271 0 : }
272 : else {
273 0 : AliFatal("Invalid digitization task! Exiting!");
274 : }
275 0 : }
276 :
277 : //____________________________________________________________________________
278 : void AliADDigitizer::DigitizeHits()
279 : {
280 : // Digitize the hits to the level of
281 : // SDigits (fTime arrays)
282 0 : Int_t nTotPhot[16];
283 0 : Float_t PMTime[16];
284 0 : Float_t PMTimeWeight[16];
285 0 : Int_t nPMHits[16];
286 :
287 0 : for(Int_t i = 0 ; i < 16; ++i) {
288 0 : memset(fTime[i],0,fNBins[i]*sizeof(Float_t));
289 0 : fLabels[i][0] = fLabels[i][1] = fLabels[i][2] = -1;
290 0 : nTotPhot[i] = 0;
291 0 : PMTime[i] = 10000;
292 0 : PMTimeWeight[i] = 0;
293 0 : nPMHits[i] = 0;
294 : }
295 :
296 0 : AliLoader* loader = fAD->GetLoader();
297 0 : if (!loader) {
298 0 : AliError("Can not get AD Loader!");
299 0 : return;
300 : }
301 0 : loader->LoadHits();
302 0 : TTree* treeH = loader->TreeH();
303 0 : if (!treeH) {
304 0 : AliError("Cannot get TreeH!");
305 0 : return;
306 : }
307 0 : TClonesArray* hits = fAD->Hits();
308 :
309 : //Loop over hits
310 0 : Int_t nTracks = (Int_t) treeH->GetEntries();
311 0 : for(Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
312 0 : fAD->ResetHits();
313 0 : treeH->GetEvent(iTrack);
314 0 : Int_t nHits = hits->GetEntriesFast();
315 0 : for (Int_t iHit = 0; iHit < nHits; iHit++) {
316 0 : AliADhit* hit = (AliADhit *)hits->UncheckedAt(iHit);
317 0 : Int_t nPhot = hit->GetNphot();
318 0 : Int_t pmt = hit->GetCell();
319 0 : if (pmt < 0) continue;
320 : //Cut late hits
321 0 : if(pmt<8 && TMath::Abs(hit->GetTof()-65.19)>3)continue;
322 0 : if(pmt>7 && TMath::Abs(hit->GetTof()-56.7 )>3)continue;
323 0 : Int_t trackLabel = hit->GetTrack();
324 0 : for(Int_t l = 0; l < 3; ++l) {
325 0 : if (fLabels[pmt][l] < 0) {
326 0 : fLabels[pmt][l] = trackLabel;
327 0 : break;
328 : }
329 : }
330 0 : Float_t dt_scintillator = gRandom->Gaus(0,kADIntTimeRes);
331 0 : Float_t t = dt_scintillator + hit->GetTof();
332 0 : nTotPhot[pmt] += nPhot;
333 0 : nPMHits[pmt]++;
334 : //PMTime[pmt] += t*nPhot*nPhot;
335 : //PMTimeWeight[pmt] += nPhot*nPhot;
336 0 : if(PMTime[pmt]>t)PMTime[pmt] = t;
337 :
338 0 : }//hit loop
339 : }//track loop
340 :
341 : //Now makes SDigits from hits
342 0 : for(Int_t iPM = 0; iPM < 16; iPM++) {
343 0 : if(nPMHits[iPM]==0 || nTotPhot[iPM]==0){
344 0 : PMTime[iPM] = 0.0;
345 0 : continue;
346 : }
347 : //PMTime[iPM] = PMTime[iPM]/PMTimeWeight[iPM];
348 0 : PMTime[iPM] += fHptdcOffset[iPM];
349 :
350 0 : fChargeSignalShape->SetParameters(fCssOffset[iPM],fCssTau[iPM],fCssSigma[iPM]);
351 0 : Float_t integral = fChargeSignalShape->Integral(0,300);
352 : //std::cout<<"Integral = "<<integral<<std::endl;
353 :
354 0 : Float_t charge = nTotPhot[iPM]*fPmGain[iPM]*fBinSize[iPM]/integral;
355 :
356 0 : Int_t firstBin = TMath::Max(0,(Int_t)((PMTime[iPM])/fBinSize[iPM]));
357 0 : Int_t lastBin = fNBins[iPM]-1;
358 : //std::cout<<"First Bin: "<<firstBin<<std::endl;
359 0 : for(Int_t iBin = firstBin; iBin <= lastBin; ++iBin) {
360 0 : Float_t tempT = fBinSize[iPM]*(0.5+iBin)-PMTime[iPM];
361 0 : if(tempT<=0)continue;
362 0 : fTime[iPM][iBin] += charge*fChargeSignalShape->Eval(tempT);
363 0 : }
364 0 : }//PM loop
365 0 : loader->UnloadHits();
366 0 : }
367 :
368 : //____________________________________________________________________________
369 : void AliADDigitizer::DigitizeSDigits()
370 : {
371 : // Digitize the fTime arrays (SDigits) to the level of
372 : // Digits (fAdc arrays)
373 0 : Float_t fMCTime[16];
374 0 : for(Int_t i = 0 ; i < 16; ++i) {
375 0 : for(Int_t j = 0; j < kADNClocks; ++j) fAdc[i][j] = 0;
376 0 : fMCTime[i] = fLeadingTime[i] = fTimeWidth[i] = 0;
377 : }
378 :
379 :
380 0 : for (Int_t ipmt = 0; ipmt < 16; ++ipmt) {
381 :
382 0 : fChargeSignalShape->SetParameters(fCssOffset[ipmt],fCssTau[ipmt],fCssSigma[ipmt]);
383 0 : Float_t maximum = 0.9*fChargeSignalShape->GetMaximum(0,300);
384 0 : Float_t integral = fChargeSignalShape->Integral(0,300);
385 0 : Float_t thr = fCalibData->GetCalibDiscriThr(ipmt)*kADChargePerADC*maximum*fBinSize[ipmt]/integral;
386 : //Float_t thr = 0;
387 :
388 : Bool_t ltFound = kFALSE, ttFound = kFALSE;
389 0 : for (Int_t iBin = 1; iBin < fNBins[ipmt]; ++iBin) {
390 0 : Float_t t = fBinSize[ipmt]*Float_t(iBin);
391 0 : if (fTime[ipmt][iBin] > 0.0) {
392 0 : if (!ltFound && (iBin < fNBinsLT[ipmt])) {
393 : ltFound = kTRUE;
394 0 : fMCTime[ipmt] = t;
395 : //std::cout<<"Leading Bin: "<<iBin<<std::endl;
396 : //std::cout<<"Leading TADC: "<<t-fClockOffset[ipmt]<<std::endl;
397 0 : }
398 : }
399 0 : if(fTime[ipmt][iBin-1] > thr && fTime[ipmt][iBin] < thr){
400 0 : if (ltFound) {
401 0 : if (!ttFound) {
402 : ttFound = kTRUE;
403 0 : fTimeWidth[ipmt] = t - fMCTime[ipmt];
404 0 : }
405 : }
406 : }
407 0 : Float_t tadc = t - fClockOffset[ipmt];
408 0 : Int_t clock = kADNClocks/2 + Int_t(tadc/25.0);
409 0 : if (clock >= 0 && clock < kADNClocks)
410 0 : fAdc[ipmt][clock] += fTime[ipmt][iBin]/kADChargePerADC;
411 : }
412 0 : AliDebug(1,Form("Channel %d Offset %f Time %f",ipmt,fClockOffset[ipmt],fMCTime[ipmt]));
413 0 : Int_t board = AliADCalibData::GetBoardNumber(ipmt);
414 0 : if (ltFound && ttFound) {
415 0 : fTimeWidth[ipmt] = fCalibData->GetWidthResolution(board)*
416 0 : Float_t(Int_t(fTimeWidth[ipmt]/fCalibData->GetWidthResolution(board)));
417 0 : if (fTimeWidth[ipmt] < Float_t(kADMinTDCWidth)*fCalibData->GetWidthResolution(board))
418 0 : fTimeWidth[ipmt] = Float_t(kADMinTDCWidth)*fCalibData->GetWidthResolution(board);
419 0 : if (fTimeWidth[ipmt] > Float_t(kADMaxTDCWidth)*fCalibData->GetWidthResolution(board))
420 0 : fTimeWidth[ipmt] = Float_t(kADMaxTDCWidth)*fCalibData->GetWidthResolution(board);
421 : }
422 : }
423 :
424 0 : fEvenOrOdd = gRandom->Integer(2);
425 0 : for (Int_t j=0; j<16; ++j){
426 : Float_t adcSignal = 0.0;
427 : Float_t adcClock = 0.0;
428 0 : for (Int_t iClock = 0; iClock < kADNClocks; ++iClock) {
429 0 : Int_t integrator = (iClock + fEvenOrOdd) % 2;
430 0 : AliDebug(1,Form("ADC %d %d %f",j,iClock,fAdc[j][iClock]));
431 0 : fAdc[j][iClock] += gRandom->Gaus(fAdcPedestal[j][integrator], fAdcSigma[j][integrator]);
432 : }
433 0 : for (Int_t iClock = 0; iClock < kADNClocks; ++iClock) {
434 0 : Int_t integrator = (iClock + fEvenOrOdd) % 2;
435 0 : adcClock = (Int_t)fAdc[j][iClock];
436 0 : if(fAdc[j][iClock]>1023) adcClock = 1023;
437 0 : adcClock -= fAdcPedestal[j][integrator];
438 0 : if(adcClock< 4*fAdcSigma[j][integrator]) adcClock = 0;
439 0 : adcSignal += adcClock;
440 : }
441 0 : fThresholdShape->SetParameter(0,fCalibData->GetCalibDiscriThr(j));
442 0 : if(gRandom->Rndm() > fThresholdShape->Eval(adcSignal)) fMCTime[j] = -1024.0;
443 0 : if(fThresholdShape->Eval(adcSignal)<1e-2) fMCTime[j] = -1024.0;
444 0 : fLeadingTime[j] = UnCorrectLeadingTime(j,fMCTime[j],adcSignal);
445 :
446 : }
447 : //Fill BB and BG flags in trigger simulator
448 0 : AliADTriggerSimulator * triggerSimulator = new AliADTriggerSimulator();
449 0 : triggerSimulator->FillFlags(fBBFlag,fBGFlag,fLeadingTime);
450 :
451 0 : }
452 :
453 : //____________________________________________________________________________
454 : void AliADDigitizer::ReadSDigits()
455 : {
456 : // Read SDigits which are then to precessed
457 : // in the following method
458 0 : for(Int_t i = 0 ; i < 16; ++i) {
459 0 : memset(fTime[i],0,fNBins[i]*sizeof(Float_t));
460 0 : fLabels[i][0] = fLabels[i][1] = fLabels[i][2] = -1;
461 : }
462 :
463 : // Loop over input files
464 0 : Int_t nFiles= fDigInput->GetNinputs();
465 0 : for (Int_t inputFile = 0; inputFile < nFiles; inputFile++) {
466 : // Get the current loader
467 : AliRunLoader* currentLoader =
468 0 : AliRunLoader::GetRunLoader(fDigInput->GetInputFolderName(inputFile));
469 :
470 0 : AliLoader *loader = currentLoader->GetLoader("ADLoader");
471 0 : loader->LoadSDigits("READ");
472 :
473 : // Get the tree of summable digits
474 0 : TTree* sdigitsTree = loader->TreeS();
475 0 : if (!sdigitsTree) {
476 0 : AliError("No sdigit tree from digInput");
477 0 : continue;
478 : }
479 :
480 : // Get the branch
481 0 : TBranch* sdigitsBranch = sdigitsTree->GetBranch("ADSDigit");
482 0 : if (!sdigitsBranch) {
483 0 : AliError("Failed to get sdigit branch");
484 0 : return;
485 : }
486 :
487 : // Set the branch address
488 0 : TClonesArray *sdigitsArray = NULL;
489 0 : sdigitsBranch->SetAddress(&sdigitsArray);
490 :
491 : // Sum contributions from the sdigits
492 : // Get number of entries in the tree
493 0 : Int_t nentries = Int_t(sdigitsBranch->GetEntries());
494 0 : for (Int_t entry = 0; entry < nentries; ++entry) {
495 0 : sdigitsBranch->GetEntry(entry);
496 : // Get the number of sdigits
497 0 : Int_t nsdigits = sdigitsArray->GetEntries();
498 :
499 0 : for (Int_t sdigit = 0; sdigit < nsdigits; sdigit++) {
500 0 : AliADSDigit* sDigit = static_cast<AliADSDigit*>(sdigitsArray->UncheckedAt(sdigit));
501 0 : Int_t pmNumber = sDigit->PMNumber();
502 0 : Int_t nbins = sDigit->GetNBins();
503 0 : if (nbins != fNBins[pmNumber]) {
504 0 : AliError(Form("Incompatible number of bins between digitizer (%d) and sdigit (%d) for PM %d! Skipping sdigit!",
505 : fNBins[pmNumber],nbins,pmNumber));
506 0 : continue;
507 : }
508 : // Sum the charges
509 0 : Float_t *charges = sDigit->GetCharges();
510 0 : for(Int_t iBin = 0; iBin < nbins; ++iBin) fTime[pmNumber][iBin] += charges[iBin];
511 : // and the labels
512 0 : Int_t *labels = sDigit->GetTracks();
513 : Int_t j = 0;
514 0 : for(Int_t i = 0; i < 3; ++i) {
515 0 : if (fLabels[pmNumber][i] < 0) {
516 0 : if (labels[j] < 0) break;
517 0 : fLabels[pmNumber][i] = labels[j];
518 0 : j++;
519 0 : }
520 : }
521 0 : }
522 : }
523 0 : loader->UnloadSDigits();
524 0 : }
525 0 : }
526 :
527 :
528 : //____________________________________________________________________________
529 : void AliADDigitizer::WriteDigits(AliLoader *loader)
530 : {
531 : // Take fAdc arrays filled by the previous
532 : // method and produce and add digits to the digit Tree
533 :
534 0 : loader->LoadDigits("UPDATE");
535 :
536 0 : if (!loader->TreeD()) loader->MakeTree("D");
537 0 : loader->MakeDigitsContainer();
538 0 : TTree* treeD = loader->TreeD();
539 0 : DigitsArray();
540 0 : treeD->Branch("ADDigit", &fDigits);
541 :
542 0 : Short_t *chargeADC = new Short_t[kADNClocks];
543 0 : for (Int_t i=0; i<16; i++) {
544 0 : for (Int_t j = 0; j < kADNClocks; ++j) {
545 0 : Int_t tempadc = Int_t(fAdc[i][j]);
546 0 : if (tempadc > 1023) tempadc = 1023;
547 0 : chargeADC[j] = tempadc;
548 : }
549 0 : AddDigit(i, fLeadingTime[i], fTimeWidth[i], Bool_t((10+fEvenOrOdd)%2), chargeADC, fBBFlag[i], fBGFlag[i], fLabels[i]);
550 : }
551 0 : delete [] chargeADC;
552 :
553 0 : treeD->Fill();
554 0 : loader->WriteDigits("OVERWRITE");
555 0 : loader->UnloadDigits();
556 0 : ResetDigits();
557 0 : }
558 :
559 : //____________________________________________________________________________
560 : void AliADDigitizer::WriteSDigits(AliLoader *loader)
561 : {
562 : // Take fTime arrays filled by the previous
563 : // method and produce and add sdigits to the sdigit Tree
564 :
565 0 : loader->LoadSDigits("UPDATE");
566 :
567 0 : if (!loader->TreeS()) loader->MakeTree("S");
568 0 : loader->MakeSDigitsContainer();
569 0 : TTree* treeS = loader->TreeS();
570 0 : SDigitsArray();
571 0 : treeS->Branch("ADSDigit", &fDigits);
572 : //fAD->MakeBranchInTree(treeS,"AD",&fDigits,8000,"");
573 :
574 0 : for (Int_t ipmt = 0; ipmt < 16; ++ipmt) {
575 0 : AddSDigit(ipmt,fNBins[ipmt],fTime[ipmt],fLabels[ipmt]);
576 : }
577 :
578 0 : treeS->Fill();
579 0 : loader->WriteSDigits("OVERWRITE");
580 0 : loader->UnloadSDigits();
581 0 : ResetDigits();
582 0 : }
583 :
584 :
585 :
586 : //____________________________________________________________________________
587 : void AliADDigitizer::AddDigit(Int_t pmnumber, Float_t time, Float_t width, Bool_t integrator, Short_t *chargeADC, Bool_t bbFlag, Bool_t bgFlag, Int_t *labels)
588 : {
589 :
590 : // Adds Digit
591 :
592 0 : TClonesArray &ldigits = *fDigits;
593 :
594 0 : new(ldigits[fNdigits++]) AliADdigit(pmnumber,time,width,integrator,chargeADC,bbFlag,bgFlag,labels);
595 :
596 0 : }
597 : //____________________________________________________________________________
598 : void AliADDigitizer::AddSDigit(Int_t pmnumber, Int_t nbins, Float_t *charges, Int_t *labels)
599 : {
600 :
601 : // Adds SDigit
602 :
603 0 : TClonesArray &ldigits = *fDigits;
604 :
605 0 : new(ldigits[fNdigits++]) AliADSDigit(pmnumber,nbins,charges,labels);
606 :
607 0 : }
608 : //____________________________________________________________________________
609 : void AliADDigitizer::ResetDigits()
610 : {
611 :
612 : // Clears Digits
613 :
614 0 : fNdigits = 0;
615 0 : if (fDigits) fDigits->Clear();
616 0 : }
617 :
618 : //____________________________________________________________________________
619 : AliADCalibData* AliADDigitizer::GetCalibData() const
620 :
621 : {
622 0 : AliCDBManager *man = AliCDBManager::Instance();
623 :
624 : AliCDBEntry *entry=0;
625 :
626 0 : entry = man->Get("AD/Calib/Data");
627 0 : if(!entry){
628 0 : AliWarning("Load of calibration data from default storage failed!");
629 0 : AliWarning("Calibration data will be loaded from local storage ($ALICE_ROOT)");
630 :
631 0 : man->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
632 0 : man->SetRun(1);
633 0 : entry = man->Get("AD/Calib/Data");
634 0 : }
635 : // Retrieval of data in directory AD/Calib/Data:
636 :
637 : AliADCalibData *calibdata = 0;
638 :
639 0 : if (entry) calibdata = (AliADCalibData*) entry->GetObject();
640 0 : if (!calibdata) AliFatal("No calibration data from calibration database !");
641 :
642 : //calibdata->PrintConfig();
643 0 : return calibdata;
644 :
645 0 : }
646 : //____________________________________________________________________________
647 : Float_t AliADDigitizer::UnCorrectLeadingTime(Int_t i, Float_t time, Float_t adc) const
648 : {
649 : // UnCorrect the MC time
650 : // for slewing effect and
651 : // misalignment of the channels
652 : const Double_t fTOF[4] = {65.2418, 65.1417, 56.6459, 56.7459};
653 :
654 0 : if (time < 1e-6) return time;
655 0 : if (adc < 1) return time;
656 :
657 : // Slewing and offset correction
658 0 : Int_t board = AliADCalibData::GetBoardNumber(i);
659 : //std::cout<<"MC time: "<<time<<std::endl;
660 0 : time -= fHptdcOffset[i];
661 : //std::cout<<"TOF: "<<time<<std::endl;
662 0 : time -= fTOF[i/4];
663 0 : if(adc<30 && fTimeSlewingExtpol[i]) time += fTimeSlewingExtpol[i]->Eval(TMath::Log10(1/adc))*fCalibData->GetTimeResolution(board);
664 :
665 0 : else time += fTimeSlewingSpline[i]->Eval(TMath::Log10(1/adc))*fCalibData->GetTimeResolution(board);
666 : //std::cout<<"Charge: "<<adc<<std::endl;
667 : //std::cout<<"Leading time: "<<time<<std::endl;
668 :
669 0 : Float_t smearedTime = SmearLeadingTime(i,time);
670 :
671 : return smearedTime;
672 0 : }
673 : //____________________________________________________________________________
674 : Float_t AliADDigitizer::SmearLeadingTime(Int_t i, Float_t time) const
675 : {
676 :
677 0 : Int_t runNumber = AliCDBManager::Instance()->GetRun();
678 : Float_t sigmaADA = 0, sigmaADC = 0;
679 0 : if(runNumber < 225753){sigmaADA = 0.45; sigmaADC = 0.15;}
680 0 : if(runNumber > 225753 && runNumber < 226501){sigmaADA = 1.0; sigmaADC = 0.15;}
681 0 : if(runNumber > 226501){sigmaADA = 0.50; sigmaADC = 0.15;}
682 :
683 0 : if(i<8)time += gRandom->Gaus(1.25,sigmaADC);
684 0 : else time += gRandom->Gaus(1.05,sigmaADA);
685 :
686 0 : return time;
687 : }
688 : //_____________________________________________________________________________
689 : void AliADDigitizer::GetTimeSlewingSplines()
690 : {
691 :
692 0 : AliCDBManager *man = AliCDBManager::Instance();
693 :
694 : AliCDBEntry *entry=0;
695 :
696 0 : entry = man->Get("AD/Calib/TimeSlewing");
697 :
698 : TList *fListSplines = 0;
699 :
700 0 : if (entry) fListSplines = (TList*) entry->GetObject();
701 0 : if (!fListSplines) AliFatal("No time slewing correction from calibration database !");
702 :
703 0 : for(Int_t i=0; i<16; i++) fTimeSlewingSpline[i] = (TSpline3*)(fListSplines->At(i));
704 :
705 :
706 0 : }
707 : //_____________________________________________________________________________
708 : void AliADDigitizer::ExtrapolateSplines()
709 : {
710 :
711 : TH1F *hTimeVsSignal;
712 :
713 0 : for(Int_t i=0; i<16; i++){
714 0 : TCanvas *c = new TCanvas("c", " ",0,0,1,1);
715 0 : c->cd();
716 0 : fTimeSlewingSpline[i]->Paint();
717 0 : hTimeVsSignal = fTimeSlewingSpline[i]->GetHistogram();
718 :
719 0 : TString TimeSlewingFitName = "hTimeSlewingFit";
720 0 : TimeSlewingFitName += i;
721 0 : fTimeSlewingExtpol[i] = new TF1(TimeSlewingFitName.Data(),"[0]+[1]*TMath::Power(10,-x*[2])",-3,0);
722 0 : fTimeSlewingExtpol[i]->SetParameter(0,650);
723 0 : fTimeSlewingExtpol[i]->SetParLimits(0,200,3000);
724 0 : fTimeSlewingExtpol[i]->SetParameter(1,450);
725 0 : fTimeSlewingExtpol[i]->SetParLimits(1,50,1000);
726 0 : fTimeSlewingExtpol[i]->SetParameter(2,-0.5);
727 0 : fTimeSlewingExtpol[i]->SetParLimits(2,-0.9,-0.05);
728 0 : fTimeSlewingExtpol[i]->SetLineColor(kMagenta);
729 0 : Int_t fitStatus = hTimeVsSignal->Fit(TimeSlewingFitName.Data(),"R"," ",-2.5,-1.5);
730 0 : if(fitStatus != 0) {
731 0 : AliWarning(Form("Extrapolation of spline %d not succesfull",i));
732 0 : fTimeSlewingExtpol[i] = 0x0;
733 0 : }
734 0 : delete c;
735 0 : }
736 0 : }
737 : //____________________________________________________________________________
738 : double AliADDigitizer::ChargeSignalShape(double *x, double *par)
739 : {
740 : // this function simulates the charge shape
741 :
742 0 : Double_t xx = x[0];
743 :
744 0 : return TMath::Exp(-0.5*TMath::Power(TMath::Log((xx+par[0])/par[1])/par[2],2));
745 : }
746 :
747 : //____________________________________________________________________________
748 : double AliADDigitizer::ThresholdShape(double *x, double *par)
749 : {
750 : // this function simulates the threshold shape
751 :
752 0 : Double_t xx = x[0];
753 :
754 0 : return 1/(1+TMath::Exp(-xx + par[0]));
755 : }
756 :
757 : //____________________________________________________________________________
758 : double AliADDigitizer::TimeSignalShape(double *x, double *par)
759 : {
760 : // this function simulates the time shape
761 :
762 0 : Double_t xx = x[0];
763 0 : if (xx <= par[0]) return 0;
764 0 : Double_t a = 1./TMath::Power((xx-par[0])/par[1],1./par[2]);
765 0 : if (xx <= par[3]) return a;
766 0 : Double_t b = 1./TMath::Power((xx-par[3])/par[4],1./par[5]);
767 0 : Double_t f = a*b/(a+b);
768 0 : AliDebug(100,Form("x=%f func=%f",xx,f));
769 : return f;
770 0 : }
771 : //____________________________________________________________________
772 : TClonesArray* AliADDigitizer::DigitsArray()
773 : {
774 : // Initialize digit array if not already and
775 : // return pointer to it.
776 0 : if (!fDigits) {
777 0 : fDigits = new TClonesArray("AliADdigit", 16);
778 0 : fNdigits = 0;
779 0 : }
780 0 : return fDigits;
781 0 : }
782 :
783 : //____________________________________________________________________
784 : TClonesArray* AliADDigitizer::SDigitsArray()
785 : {
786 : // Initialize sdigit array if not already and
787 : // return pointer to it.
788 0 : if (!fDigits) {
789 0 : fDigits = new TClonesArray("AliADSDigit", 16);
790 0 : fNdigits = 0;
791 0 : }
792 0 : return fDigits;
793 0 : }
794 :
|
