]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PHOS/AliPHOSDigitizer.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Added DCA to the THn (M Tangaro)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSDigitizer.cxx
... / ...
CommitLineData
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$ */
17
18/* History of cvs commits:
19 *
20 * $Log: AliPHOSDigitizer.cxx,v $
21 * Revision 1.104 2007/12/18 09:08:18 hristov
22 * Splitting of the QA maker into simulation and reconstruction dependent parts (Yves)
23 *
24 * Revision 1.103 2007/11/07 11:25:06 schutz
25 * Comment out the QA checking before starting digitization
26 *
27 * Revision 1.102 2007/10/19 18:04:29 schutz
28 * The standalone QA data maker is called from AliSimulation and AliReconstruction outside the event loop; i.e. re-reading the data. The QA data making in the event loop has been commented out.
29 *
30 * Revision 1.101 2007/10/14 21:08:10 schutz
31 * Introduced the checking of QA results from previous step before entering the event loop
32 *
33 * Revision 1.100 2007/10/10 09:05:10 schutz
34 * Changing name QualAss to QA
35 *
36 * Revision 1.99 2007/09/30 17:08:20 schutz
37 * Introducing the notion of QA data acquisition cycle (needed by online)
38 *
39 * Revision 1.98 2007/09/26 14:22:17 cvetan
40 * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
41 *
42 * Revision 1.97 2007/08/07 14:12:03 kharlov
43 * Quality assurance added (Yves Schutz)
44 *
45 * Revision 1.96 2007/04/28 10:43:36 policheh
46 * Dead channels simulation: digit energy sets to 0.
47 *
48 * Revision 1.95 2007/04/10 07:20:52 kharlov
49 * Decalibration should use the same CDB as calibration in AliPHOSClusterizerv1
50 *
51 * Revision 1.94 2007/02/01 10:34:47 hristov
52 * Removing warnings on Solaris x86
53 *
54 * Revision 1.93 2006/10/17 13:17:01 kharlov
55 * Replace AliInfo by AliDebug
56 *
57 * Revision 1.92 2006/08/28 10:01:56 kharlov
58 * Effective C++ warnings fixed (Timur Pocheptsov)
59 *
60 * Revision 1.91 2006/04/29 20:25:30 hristov
61 * Decalibration is implemented (Yu.Kharlov)
62 *
63 * Revision 1.90 2006/04/22 10:30:17 hristov
64 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
65 *
66 * Revision 1.89 2006/04/11 15:22:59 hristov
67 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
68 *
69 * Revision 1.88 2006/03/13 14:05:43 kharlov
70 * Calibration objects for EMC and CPV
71 *
72 * Revision 1.87 2005/08/24 15:33:49 kharlov
73 * Calibration data for raw digits
74 *
75 * Revision 1.86 2005/07/12 20:07:35 hristov
76 * Changes needed to run simulation and reconstrruction in the same AliRoot session
77 *
78 * Revision 1.85 2005/05/28 14:19:04 schutz
79 * Compilation warnings fixed by T.P.
80 *
81 */
82
83//_________________________________________________________________________
84//*-- Author : Dmitri Peressounko (SUBATECH & Kurchatov Institute)
85//////////////////////////////////////////////////////////////////////////////
86// This class performs digitization of Summable digits (in the PHOS case it is just
87// the sum of contributions from all primary particles into a given cell).
88// In addition it performs mixing of summable digits from different events.
89// The name of the class is also the title of the branch that will contain
90// the created SDigits
91// The title of the class is the name of the file that contains the hits from
92// which the SDigits are created
93//
94// For each event two branches are created in TreeD:
95// "PHOS" - list of digits
96// "AliPHOSDigitizer" - AliPHOSDigitizer with all parameters used in digitization
97//
98// Note, that one can set a title for new digits branch, and repeat digitization with
99// another set of parameters.
100//
101// Use case:
102// root[0] AliPHOSDigitizer * d = new AliPHOSDigitizer() ;
103// root[1] d->Digitize()
104// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
105// //Digitizes SDigitis in all events found in file galice.root
106//
107// root[2] AliPHOSDigitizer * d1 = new AliPHOSDigitizer("galice1.root") ;
108// // Will read sdigits from galice1.root
109// root[3] d1->MixWith("galice2.root")
110// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
111// // Reads another set of sdigits from galice2.root
112// root[3] d1->MixWith("galice3.root")
113// // Reads another set of sdigits from galice3.root
114// root[4] d->Digitize("deb timing")
115// // Reads SDigits from files galice1.root, galice2.root ....
116// // mixes them and stores produced Digits in file galice1.root
117// // deb - prints number of produced digits
118// // deb all - prints list of produced digits
119// // timing - prints time used for digitization
120//
121
122// --- ROOT system ---
123#include "TTree.h"
124#include "TSystem.h"
125#include "TBenchmark.h"
126#include "TRandom.h"
127#include "TMath.h"
128
129// --- Standard library ---
130
131// --- AliRoot header files ---
132#include <TGeoManager.h>
133#include "AliLog.h"
134#include "AliDigitizationInput.h"
135#include "AliPHOSDigit.h"
136#include "AliPHOSDigitizer.h"
137#include "AliPHOSGeometry.h"
138#include "AliPHOSTick.h"
139#include "AliPHOSSimParam.h"
140#include "AliPHOSCalibData.h"
141#include "AliRunLoader.h"
142#include "AliPHOSLoader.h"
143#include "AliPHOSPulseGenerator.h"
144
145ClassImp(AliPHOSDigitizer)
146
147
148//____________________________________________________________________________
149AliPHOSDigitizer::AliPHOSDigitizer() :
150 AliDigitizer("",""),
151 fDefaultInit(kTRUE),
152 fDigitsInRun(0),
153 fInit(kFALSE),
154 fInput(0),
155 fInputFileNames(0x0),
156 fEventNames(0x0),
157 fEmcCrystals(0),
158 fEventFolderName(""),
159 fFirstEvent(0),
160 fLastEvent(0),
161 fcdb(0x0),
162 fEventCounter(0),
163 fPulse(0),
164 fADCValuesLG(0),
165 fADCValuesHG(0)
166{
167 // ctor
168 InitParameters() ;
169 fDigInput = 0 ; // We work in the standalong mode
170}
171
172//____________________________________________________________________________
173AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName,
174 TString eventFolderName):
175 AliDigitizer("PHOSDigitizer", alirunFileName),
176 fDefaultInit(kFALSE),
177 fDigitsInRun(0),
178 fInit(kFALSE),
179 fInput(0),
180 fInputFileNames(0x0),
181 fEventNames(0x0),
182 fEmcCrystals(0),
183 fEventFolderName(eventFolderName),
184 fFirstEvent(0),
185 fLastEvent(0),
186 fcdb(0x0),
187 fEventCounter(0),
188 fPulse(0),
189 fADCValuesLG(0),
190 fADCValuesHG(0)
191{
192 // ctor
193 InitParameters() ;
194 Init() ;
195 fDefaultInit = kFALSE ;
196 fDigInput = 0 ; // We work in the standalone mode
197 fcdb = new AliPHOSCalibData(-1);
198}
199
200//____________________________________________________________________________
201AliPHOSDigitizer::AliPHOSDigitizer(AliDigitizationInput * rd) :
202 AliDigitizer(rd,"PHOSDigitizer"),
203 fDefaultInit(kFALSE),
204 fDigitsInRun(0),
205 fInit(kFALSE),
206 fInput(0),
207 fInputFileNames(0x0),
208 fEventNames(0x0),
209 fEmcCrystals(0),
210 fEventFolderName(fDigInput->GetInputFolderName(0)),
211 fFirstEvent(0),
212 fLastEvent(0),
213 fcdb (0x0),
214 fEventCounter(0),
215 fPulse(0),
216 fADCValuesLG(0),
217 fADCValuesHG(0)
218
219{
220 // ctor Init() is called by RunDigitizer
221 fDigInput = rd ;
222 SetTitle(static_cast<AliStream*>(fDigInput->GetInputStream(0))->GetFileName(0));
223 InitParameters() ;
224 fDefaultInit = kFALSE ;
225 fcdb = new AliPHOSCalibData(-1);
226}
227
228//____________________________________________________________________________
229 AliPHOSDigitizer::~AliPHOSDigitizer()
230{
231 // dtor
232 delete [] fInputFileNames ;
233 delete [] fEventNames ;
234
235 delete fPulse;
236 delete [] fADCValuesLG;
237 delete [] fADCValuesHG;
238
239 if(fcdb){ delete fcdb ; fcdb=0;}
240
241}
242
243//____________________________________________________________________________
244void AliPHOSDigitizer::Digitize(Int_t event)
245{
246
247 // Makes the digitization of the collected summable digits.
248 // It first creates the array of all PHOS modules
249 // filled with noise (different for EMC, and CPV) and
250 // then adds contributions from SDigits.
251 // This design avoids scanning over the list of digits to add
252 // contribution to new SDigits only.
253
254 Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
255
256 //First stream
257 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
258 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
259
260 Int_t readEvent = event ;
261 if (fDigInput)
262 readEvent = static_cast<AliStream*>(fDigInput->GetInputStream(0))->GetCurrentEventNumber() ;
263 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
264 readEvent, GetTitle(), fEventFolderName.Data())) ;
265 rl->GetEvent(readEvent) ;
266 phosLoader->CleanSDigits() ;
267 phosLoader->LoadSDigits("READ") ;
268
269 //Prepare Output
270 TClonesArray * digits = phosLoader->Digits() ;
271 if( !digits ) {
272 phosLoader->MakeDigitsArray() ;
273 digits = phosLoader->Digits() ;
274 }
275
276 digits->Clear() ;
277
278 //
279 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
280 //Making digits with noise, first EMC
281 //Check which PHOS modules are present
282 Bool_t isPresent[5] ;
283 TString volpath ;
284 Int_t nmod=0 ;
285 for(Int_t i=0; i<5; i++){
286 volpath = "/ALIC_1/PHOS_";
287 volpath += i+1;
288 if (gGeoManager->CheckPath(volpath.Data())) {
289 isPresent[i]=1 ;
290 nmod++ ;
291 }
292 else{
293 isPresent[i]=0 ;
294 }
295 }
296
297 Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
298
299 Int_t nCPV ;
300 Int_t absID ;
301
302 //check if CPV exists
303 Bool_t isCPVpresent=0 ;
304 for(Int_t i=1; i<=5 && !isCPVpresent; i++){
305 volpath = "/ALIC_1/PHOS_";
306 volpath += i;
307 volpath += "/PCPV_1";
308 if (gGeoManager->CheckPath(volpath.Data()))
309 isCPVpresent=1 ;
310 }
311
312 if(isCPVpresent){
313 nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
314 }
315 else{
316 nCPV = nEMC ;
317 }
318
319 digits->Expand(nCPV) ;
320
321 //take all the inputs to add together and load the SDigits
322 TObjArray * sdigArray = new TObjArray(fInput) ;
323 sdigArray->AddAt(phosLoader->SDigits(), 0) ;
324
325 for(Int_t i = 1 ; i < fInput ; i++){
326 TString tempo(fEventNames[i]) ;
327 tempo += i ;
328 AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
329 if(!rl2){
330 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
331 if (!rl2) {
332 Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
333 return ;
334 }
335 rl2->LoadHeader();
336 }
337 AliPHOSLoader * phosLoader2 = static_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
338
339 if(fDigInput){
340 readEvent = static_cast<AliStream*>(fDigInput->GetInputStream(i))->GetCurrentEventNumber() ;
341 }
342 TClonesArray * digs ;
343 if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
344 AliInfo(Form("Adding event %d from input stream %d %s %s",
345 readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
346 rl2->GetEvent(readEvent) ;
347 phosLoader2->CleanDigits() ;
348 phosLoader2->LoadDigits("READ") ;
349 digs = phosLoader2->Digits() ;
350 toMakeNoise=0 ; //Do not add noise, it is already in stream
351 }
352 else{
353 AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
354 readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
355 rl2->GetEvent(readEvent) ;
356 phosLoader2->CleanSDigits() ;
357 phosLoader2->LoadSDigits("READ") ;
358 digs = phosLoader2->SDigits() ;
359 }
360 sdigArray->AddAt(digs, i) ;
361 }
362
363 //Find the first crystal with signal
364 Int_t nextSig = 200000 ;
365 TClonesArray * sdigits ;
366 for(Int_t i = 0 ; i < fInput ; i++){
367 sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
368 if ( !sdigits->GetEntriesFast() )
369 continue ;
370 Int_t curNext = static_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
371 if(curNext < nextSig)
372 nextSig = curNext ;
373 }
374
375 TArrayI index(fInput) ;
376 index.Reset() ; //Set all indexes to zero
377
378 AliPHOSDigit * digit ;
379 AliPHOSDigit * curSDigit ;
380
381// TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
382
383 //Put Noise contribution
384 Double_t apdNoise = 0. ;
385 if(toMakeNoise)
386 apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ;
387
388 Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
389 Int_t idigit= 0;
390 for(Int_t imod=0; imod<5; imod++){
391 if(!isPresent[imod])
392 continue ;
393 Int_t firstAbsId=imod*emcpermod+1 ;
394 Int_t lastAbsId =(imod+1)*emcpermod ;
395 for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
396 Float_t noise = gRandom->Gaus(0.,apdNoise) ;
397 new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
398 //look if we have to add signal?
399 digit = static_cast<AliPHOSDigit *>(digits->At(idigit)) ;
400 idigit++ ;
401
402 if(absID==nextSig){
403 //Add SDigits from all inputs
404// ticks->Clear() ;
405// Int_t contrib = 0 ;
406
407//New Timing model is necessary
408// Float_t a = digit->GetEnergy() ;
409// Float_t b = TMath::Abs( a / fTimeSignalLength) ;
410// //Mark the beginning of the signal
411// new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
412// //Mark the end of the signal
413// new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
414
415// Calculate time as time of the largest digit
416 Float_t time = digit->GetTime() ;
417 Float_t eTime= digit->GetEnergy() ;
418
419 //loop over inputs
420 for(Int_t i = 0 ; i < fInput ; i++){
421 if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
422 curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
423 if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
424 curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
425 curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
426 }
427 }
428 else
429 curSDigit = 0 ;
430 //May be several digits will contribute from the same input
431 while(curSDigit && curSDigit->GetId() == absID){
432 //Shift primary to separate primaries belonging different inputs
433 Int_t primaryoffset ;
434 if(fDigInput)
435 primaryoffset = fDigInput->GetMask(i) ;
436 else
437 primaryoffset = 10000000*i ;
438 curSDigit->ShiftPrimary(primaryoffset) ;
439
440//New Timing model is necessary
441// a = curSDigit->GetEnergy() ;
442// b = a /fTimeSignalLength ;
443// new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
444// new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
445 if(curSDigit->GetEnergy()>eTime){
446 eTime=curSDigit->GetEnergy() ;
447 time=curSDigit->GetTime() ;
448 }
449 *digit += *curSDigit ; //add energies
450
451 index[i]++ ;
452 if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
453 curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
454 else
455 curSDigit = 0 ;
456 }
457 }
458
459 //calculate and set time
460//New Timing model is necessary
461// Float_t time = FrontEdgeTime(ticks) ;
462 digit->SetTime(time) ;
463
464 //Find next signal module
465 nextSig = 200000 ;
466 for(Int_t i = 0 ; i < fInput ; i++){
467 sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
468 Int_t curNext = nextSig ;
469 if(sdigits->GetEntriesFast() > index[i] ){
470 curNext = static_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
471 }
472 if(curNext < nextSig) nextSig = curNext ;
473 }
474 }
475 }
476 }
477
478
479 //Apply non-linearity
480 if(AliPHOSSimParam::GetInstance()->IsCellNonlinearityOn()){ //Apply non-lineairyt on cell level
481 const Double_t aNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyA() ;
482 const Double_t bNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyB() ;
483 const Double_t cNL = AliPHOSSimParam::GetInstance()->GetCellNonLineairyC() ;
484 for(Int_t i = 0 ; i < nEMC ; i++){
485 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
486 Double_t e= digit->GetEnergy() ;
487 // version(1) digit->SetEnergy(e*(1+a*TMath::Exp(-e/b))) ;
488 digit->SetEnergy(e*cNL*(1.+aNL*TMath::Exp(-e*e/2./bNL/bNL))) ; //Better agreement with data...
489 }
490 }
491
492
493 //distretize energy if necessary
494 if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
495 Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
496 for(Int_t i = 0 ; i < nEMC ; i++){
497 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
498 digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
499 }
500 }
501
502 //Apply decalibration if necessary
503 for(Int_t i = 0 ; i < nEMC ; i++){
504 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
505 Decalibrate(digit) ;
506 }
507
508// ticks->Delete() ;
509// delete ticks ;
510
511 //Now CPV digits (different noise and no timing)
512 Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
513 Int_t nEMCtotal=emcpermod*5 ;
514 Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
515 if(isCPVpresent){ //CPV is present in current geometry
516 for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
517 if(!isPresent[imod])
518 continue ;
519 Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
520 Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
521 for(absID = firstAbsId; absID <= lastAbsId; absID++){
522 Float_t noise = gRandom->Gaus(0., cpvNoise) ;
523 new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
524 idigit++ ;
525 //look if we have to add signal?
526 if(absID==nextSig){
527 digit = static_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
528 //Add SDigits from all inputs
529 for(Int_t i = 0 ; i < fInput ; i++){
530 if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
531 curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
532 else
533 curSDigit = 0 ;
534
535 //May be several digits will contribute from the same input
536 while(curSDigit && curSDigit->GetId() == absID){
537 //Shift primary to separate primaries belonging different inputs
538 Int_t primaryoffset ;
539 if(fDigInput)
540 primaryoffset = fDigInput->GetMask(i) ;
541 else
542 primaryoffset = 10000000*i ;
543 curSDigit->ShiftPrimary(primaryoffset) ;
544
545 //add energies
546 *digit += *curSDigit ;
547 index[i]++ ;
548 if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
549 curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
550 else
551 curSDigit = 0 ;
552 }
553 }
554
555 //Find next signal module
556 nextSig = 200000 ;
557 for(Int_t i = 0 ; i < fInput ; i++){
558 sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
559 Int_t curNext = nextSig ;
560 if(sdigits->GetEntriesFast() > index[i] )
561 curNext = static_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
562 if(curNext < nextSig) nextSig = curNext ;
563 }
564
565 }
566 }
567 }
568 }
569
570 delete sdigArray ; //We should not delete its contents
571
572 Int_t relId[4];
573
574 //set amplitudes in bad channels to zero
575
576 for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
577 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
578 geom->AbsToRelNumbering(digit->GetId(),relId);
579 if(relId[1] == 0) // Emc
580 if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.);
581 }
582
583 //remove digits below thresholds
584 Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
585 for(Int_t i = 0 ; i < nEMC ; i++){
586 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
587
588 if(digit->GetEnergy() < emcThreshold){
589 digits->RemoveAt(i) ;
590 continue ;
591 }
592
593 geom->AbsToRelNumbering(digit->GetId(),relId);
594
595// digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
596
597 Float_t tres = TimeResolution(digit->GetEnergy()) ;
598 digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
599
600 fPulse->Reset();
601 fPulse->SetAmplitude(digit->GetEnergy()/
602 fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
603 fPulse->SetTZero(digit->GetTimeR());
604 fPulse->MakeSamples();
605 fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ;
606 Int_t nSamples = fPulse->GetRawFormatTimeBins();
607 digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
608 digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
609 }
610
611 Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
612 for(Int_t i = nEMC; i < nCPV ; i++){
613 if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
614 digits->RemoveAt(i) ;
615 }
616
617 digits->Compress() ;
618 Int_t ndigits = digits->GetEntriesFast() ;
619
620 //Set indexes in list of digits and make true digitization of the energy
621 for (Int_t i = 0 ; i < ndigits ; i++) {
622 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
623 digit->SetIndexInList(i) ;
624 if(digit->GetId() > fEmcCrystals){ //digitize CPV only
625 digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
626 }
627 }
628
629}
630//____________________________________________________________________________
631Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
632 //Apply calibration
633 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
634
635 //Determine rel.position of the cell absolute ID
636 Int_t relId[4];
637 geom->AbsToRelNumbering(absId,relId);
638 Int_t module=relId[0];
639 Int_t row =relId[2];
640 Int_t column=relId[3];
641 if(relId[1]==0){ //This Is EMC
642 Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
643 return amp*calibration ;
644 }
645 return 0 ;
646}
647//____________________________________________________________________________
648void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
649{
650 // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
651
652 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
653
654 //Determine rel.position of the cell absolute ID
655 Int_t relId[4];
656 geom->AbsToRelNumbering(digit->GetId(),relId);
657 Int_t module=relId[0];
658 Int_t row =relId[2];
659 Int_t column=relId[3];
660 if(relId[1]==0){ //This Is EMC
661 Float_t decalib = fcdb->GetADCchannelEmcDecalib(module,column,row); // O(1)
662 Float_t calibration = fcdb->GetADCchannelEmc(module,column,row)*decalib;
663 Float_t energy = digit->GetEnergy()/calibration;
664 digit->SetEnergy(energy); //Now digit measures E in ADC counts
665 Float_t time = digit->GetTime() ;
666 time-=fcdb->GetTimeShiftEmc(module,column,row);
667 digit->SetTime(time) ;
668 }
669}
670//____________________________________________________________________________
671Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
672 //Apply time calibration
673 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
674
675 //Determine rel.position of the cell absolute ID
676 Int_t relId[4];
677 geom->AbsToRelNumbering(absId,relId);
678 Int_t module=relId[0];
679 Int_t row =relId[2];
680 Int_t column=relId[3];
681 time += fcdb->GetTimeShiftEmc(module,column,row);
682 return time ;
683}
684//____________________________________________________________________________
685Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
686{
687 // Returns digitized value of the CPV charge in a pad absId
688
689 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
690
691 //Determine rel.position of the cell absId
692 Int_t relId[4];
693 geom->AbsToRelNumbering(absId,relId);
694 Int_t module=relId[0];
695 Int_t row =relId[2];
696 Int_t column=relId[3];
697
698 Int_t channel = 0;
699
700 if(absId > fEmcCrystals){ //digitize CPV only
701
702 //reading calibration data for cell absId.
703 Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
704 Float_t adcChanelCpv = fcdb->GetADCchannelCpv( module,column,row);
705
706 channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
707 Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
708 if(channel > nMax ) channel = nMax ;
709 }
710 return channel ;
711}
712
713//____________________________________________________________________________
714void AliPHOSDigitizer::Digitize(Option_t *option)
715{
716 // Steering method to process digitization for events
717 // in the range from fFirstEvent to fLastEvent.
718 // This range is optionally set by SetEventRange().
719 // if fLastEvent=-1, then process events until the end.
720 // by default fLastEvent = fFirstEvent (process only one event)
721
722 if (!fInit) { // to prevent overwrite existing file
723 AliError(Form("Give a version name different from %s",
724 fEventFolderName.Data() )) ;
725 return ;
726 }
727
728 if (strstr(option,"print")) {
729 Print();
730 return ;
731 }
732
733 if(strstr(option,"tim"))
734 gBenchmark->Start("PHOSDigitizer");
735
736 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
737 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
738
739 if (fLastEvent == -1)
740 fLastEvent = rl->GetNumberOfEvents() - 1 ;
741 else if (fDigInput)
742 fLastEvent = fFirstEvent ;
743
744 Int_t nEvents = fLastEvent - fFirstEvent + 1;
745
746 Int_t ievent ;
747
748 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
749 fEventCounter++ ;
750
751 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
752
753 WriteDigits() ;
754
755 if(strstr(option,"deb"))
756 PrintDigits(option);
757
758 //increment the total number of Digits per run
759 fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
760 }
761
762 if(strstr(option,"tim")){
763 gBenchmark->Stop("PHOSDigitizer");
764 TString message ;
765 message = " took %f seconds for Digitizing %f seconds per event\n" ;
766 AliInfo(Form( message.Data(),
767 gBenchmark->GetCpuTime("PHOSDigitizer"),
768 gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
769 }
770}
771//____________________________________________________________________________
772Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
773 //calculate TOF resolution using beam-test resutls
774 Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
775 Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
776 return TMath::Sqrt(a*a+b*b/e) ;
777}
778
779////____________________________________________________________________________
780//Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
781//{
782// // Returns the shortest time among all time ticks
783//
784// ticks->Sort() ; //Sort in accordance with times of ticks
785// TIter it(ticks) ;
786// AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
787// Float_t time = ctick->CrossingTime(fTimeThreshold) ;
788//
789// AliPHOSTick * t ;
790// while((t=(AliPHOSTick*) it.Next())){
791// if(t->GetTime() < time) //This tick starts before crossing
792// *ctick+=*t ;
793// else
794// return time ;
795//
796// time = ctick->CrossingTime(fTimeThreshold) ;
797// }
798// return time ;
799//}
800
801//____________________________________________________________________________
802Bool_t AliPHOSDigitizer::Init()
803{
804 // Makes all memory allocations
805 fInit = kTRUE ;
806
807 AliPHOSGeometry *geom;
808 if (!(geom = AliPHOSGeometry::GetInstance()))
809 geom = AliPHOSGeometry::GetInstance("IHEP","");
810// const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
811
812 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
813
814 fFirstEvent = 0 ;
815 fLastEvent = fFirstEvent ;
816 if (fDigInput)
817 fInput = fDigInput->GetNinputs() ;
818 else
819 fInput = 1 ;
820
821 fInputFileNames = new TString[fInput] ;
822 fEventNames = new TString[fInput] ;
823 fInputFileNames[0] = GetTitle() ;
824 fEventNames[0] = fEventFolderName.Data() ;
825 Int_t index ;
826 for (index = 1 ; index < fInput ; index++) {
827 fInputFileNames[index] = static_cast<AliStream*>(fDigInput->GetInputStream(index))->GetFileName(0);
828 TString tempo = fDigInput->GetInputFolderName(index) ;
829 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fDigInput
830 }
831
832 //to prevent cleaning of this object while GetEvent is called
833 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
834 if(!rl){
835 AliRunLoader::Open(GetTitle(), fEventFolderName) ;
836 }
837 return fInit ;
838}
839
840//____________________________________________________________________________
841void AliPHOSDigitizer::InitParameters()
842{
843 // Set initial parameters Digitizer
844
845 fDigitsInRun = 0 ;
846 SetEventRange(0,-1) ;
847 fPulse = new AliPHOSPulseGenerator();
848 fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
849 fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
850
851}
852
853//__________________________________________________________________
854void AliPHOSDigitizer::Print(const Option_t *)const
855{
856 // Print Digitizer's parameters
857 AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
858 if( strcmp(fEventFolderName.Data(), "") != 0 ){
859 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
860
861 Int_t nStreams ;
862 if (fDigInput)
863 nStreams = GetNInputStreams() ;
864 else
865 nStreams = fInput ;
866
867 Int_t index = 0 ;
868 for (index = 0 ; index < nStreams ; index++) {
869 TString tempo(fEventNames[index]) ;
870 tempo += index ;
871 printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ;
872 }
873
874 // printf("\nWith following parameters:\n") ;
875 // printf(" Electronics noise in EMC (fPinNoise) = %f GeV\n", fPinNoise ) ;
876 // printf(" Threshold in EMC (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
877 // printf(" Noise in CPV (fCPVNoise) = %f aux units\n", fCPVNoise ) ;
878 // printf(" Threshold in CPV (fCPVDigitThreshold) = %f aux units\n",fCPVDigitThreshold ) ;
879 printf(" ---------------------------------------------------\n") ;
880 }
881 else
882 AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
883
884}
885
886//__________________________________________________________________
887 void AliPHOSDigitizer::PrintDigits(Option_t * option)
888{
889 // Print a table of digits
890
891
892 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
893 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
894 TClonesArray * digits = phosLoader->Digits() ;
895 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
896
897 AliInfo(Form("%d", digits->GetEntriesFast())) ;
898 printf("\nevent %d", gAlice->GetEvNumber()) ;
899 printf("\n Number of entries in Digits list %d", digits->GetEntriesFast() ) ;
900
901
902 if(strstr(option,"all")||strstr(option,"EMC")){
903 //loop over digits
904 AliPHOSDigit * digit;
905 printf("\nEMC digits (with primaries):\n") ;
906 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
907 Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
908 Int_t index ;
909 for (index = 0 ; (index < digits->GetEntriesFast()) &&
910 (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
911 digit = (AliPHOSDigit * ) digits->At(index) ;
912 if(digit->GetNprimary() == 0)
913 continue;
914// printf("%6d %8d %6.5e %4d %2d :",
915// digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ; // YVK
916 printf("%6d %.4f %6.5e %4d %2d :",
917 digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
918 Int_t iprimary;
919 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
920 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
921 }
922 printf("\n") ;
923 }
924 }
925
926 if(strstr(option,"all")||strstr(option,"CPV")){
927
928 //loop over CPV digits
929 AliPHOSDigit * digit;
930 printf("\nCPV digits (with primaries):\n") ;
931 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
932 Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
933 Int_t index ;
934 for (index = 0 ; index < digits->GetEntriesFast(); index++) {
935 digit = (AliPHOSDigit * ) digits->At(index) ;
936 if(digit->GetId() > maxEmc){
937 printf("%6d %8d %4d %2d :",
938 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
939 Int_t iprimary;
940 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
941 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
942 }
943 printf("\n") ;
944 }
945 }
946 }
947
948}
949
950//__________________________________________________________________
951Float_t AliPHOSDigitizer::TimeOfNoise(void) const
952{ // Calculates the time signal generated by noise
953 //PH Info("TimeOfNoise", "Change me") ;
954 return gRandom->Rndm() * 1.28E-5;
955}
956
957//__________________________________________________________________
958void AliPHOSDigitizer::Unload()
959{
960
961 Int_t i ;
962 for(i = 1 ; i < fInput ; i++){
963 TString tempo(fEventNames[i]) ;
964 tempo += i ;
965 AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ;
966 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
967 phosLoader->UnloadSDigits() ;
968 }
969
970 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
971 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
972 phosLoader->UnloadDigits() ;
973}
974
975//____________________________________________________________________________
976void AliPHOSDigitizer::WriteDigits()
977{
978
979 // Makes TreeD in the output file.
980 // Check if branch already exists:
981 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
982 // already existing branches.
983 // else creates branch with Digits, named "PHOS", title "...",
984 // and branch "AliPHOSDigitizer", with the same title to keep all the parameters
985 // and names of files, from which digits are made.
986
987 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
988 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
989
990 const TClonesArray * digits = phosLoader->Digits() ;
991 TTree * treeD = phosLoader->TreeD();
992 if(!treeD){
993 phosLoader->MakeTree("D");
994 treeD = phosLoader->TreeD();
995 }
996
997 // -- create Digits branch
998 Int_t bufferSize = 32000 ;
999 TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1000 digitsBranch->SetTitle(fEventFolderName);
1001 digitsBranch->Fill() ;
1002
1003 phosLoader->WriteDigits("OVERWRITE");
1004
1005 Unload() ;
1006
1007}