]> 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
Coverity fixes
[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 TTask 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 TTask is also the title of the branch that will contain
90// the created SDigits
91// The title of the TTAsk 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->ExecuteTask()
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->ExecuteTask("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 "AliRunDigitizer.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 fManager = 0 ; // We work in the standalong mode
170}
171
172//____________________________________________________________________________
173AliPHOSDigitizer::AliPHOSDigitizer(TString alirunFileName,
174 TString eventFolderName):
175 AliDigitizer("PHOS"+AliConfig::Instance()->GetDigitizerTaskName(), 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 fManager = 0 ; // We work in the standalone mode
197 fcdb = new AliPHOSCalibData(-1);
198}
199
200//____________________________________________________________________________
201AliPHOSDigitizer::AliPHOSDigitizer(const AliPHOSDigitizer & d) :
202 AliDigitizer(d),
203 fDefaultInit(d.fDefaultInit),
204 fDigitsInRun(d.fDigitsInRun),
205 fInit(d.fInit),
206 fInput(d.fInput),
207 fInputFileNames(0x0),//?
208 fEventNames(0x0),//?
209 fEmcCrystals(d.fEmcCrystals),
210 fEventFolderName(d.fEventFolderName),
211 fFirstEvent(d.fFirstEvent),
212 fLastEvent(d.fLastEvent),
213 fcdb (0x0),
214 fEventCounter(0),
215 fPulse(0),
216 fADCValuesLG(0),
217 fADCValuesHG(0)
218{
219 // copyy ctor
220 SetName(d.GetName()) ;
221 SetTitle(d.GetTitle()) ;
222 fcdb = new AliPHOSCalibData(-1);
223}
224
225//____________________________________________________________________________
226AliPHOSDigitizer::AliPHOSDigitizer(AliRunDigitizer * rd) :
227 AliDigitizer(rd,"PHOS"+AliConfig::Instance()->GetDigitizerTaskName()),
228 fDefaultInit(kFALSE),
229 fDigitsInRun(0),
230 fInit(kFALSE),
231 fInput(0),
232 fInputFileNames(0x0),
233 fEventNames(0x0),
234 fEmcCrystals(0),
235 fEventFolderName(fManager->GetInputFolderName(0)),
236 fFirstEvent(0),
237 fLastEvent(0),
238 fcdb (0x0),
239 fEventCounter(0),
240 fPulse(0),
241 fADCValuesLG(0),
242 fADCValuesHG(0)
243
244{
245 // ctor Init() is called by RunDigitizer
246 fManager = rd ;
247 SetTitle(static_cast<AliStream*>(fManager->GetInputStream(0))->GetFileName(0));
248 InitParameters() ;
249 fDefaultInit = kFALSE ;
250 fcdb = new AliPHOSCalibData(-1);
251}
252
253//____________________________________________________________________________
254 AliPHOSDigitizer::~AliPHOSDigitizer()
255{
256 // dtor
257 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
258 if(rl){
259 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
260
261 if(phosLoader)
262 phosLoader->CleanDigitizer() ;
263 }
264
265 delete [] fInputFileNames ;
266 delete [] fEventNames ;
267
268 delete fPulse;
269 delete [] fADCValuesLG;
270 delete [] fADCValuesHG;
271
272 if(fcdb){ delete fcdb ; fcdb=0;}
273
274}
275
276//____________________________________________________________________________
277void AliPHOSDigitizer::Digitize(Int_t event)
278{
279
280 // Makes the digitization of the collected summable digits.
281 // It first creates the array of all PHOS modules
282 // filled with noise (different for EMC, and CPV) and
283 // then adds contributions from SDigits.
284 // This design avoids scanning over the list of digits to add
285 // contribution to new SDigits only.
286
287 Bool_t toMakeNoise = kTRUE ; //Do not create noisy digits if merge with real data
288
289 //First stream
290 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
291 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
292
293 Int_t readEvent = event ;
294 if (fManager)
295 readEvent = static_cast<AliStream*>(fManager->GetInputStream(0))->GetCurrentEventNumber() ;
296 AliDebug(1,Form("Adding event %d from input stream 0 %s %s",
297 readEvent, GetTitle(), fEventFolderName.Data())) ;
298 rl->GetEvent(readEvent) ;
299 phosLoader->CleanSDigits() ;
300 phosLoader->LoadSDigits("READ") ;
301
302 //Prepare Output
303 TClonesArray * digits = phosLoader->Digits() ;
304 if( !digits ) {
305 phosLoader->MakeDigitsArray() ;
306 digits = phosLoader->Digits() ;
307 }
308
309 digits->Clear() ;
310
311 //
312 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
313 //Making digits with noise, first EMC
314 //Check which PHOS modules are present
315 Bool_t isPresent[5] ;
316 TString volpath ;
317 Int_t nmod=0 ;
318 for(Int_t i=0; i<5; i++){
319 volpath = "/ALIC_1/PHOS_";
320 volpath += i+1;
321 if (gGeoManager->CheckPath(volpath.Data())) {
322 isPresent[i]=1 ;
323 nmod++ ;
324 }
325 else{
326 isPresent[i]=0 ;
327 }
328 }
329
330 Int_t nEMC = nmod*geom->GetNPhi()*geom->GetNZ();
331
332 Int_t nCPV ;
333 Int_t absID ;
334
335 //check if CPV exists
336 Bool_t isCPVpresent=0 ;
337 for(Int_t i=1; i<=5 && !isCPVpresent; i++){
338 volpath = "/ALIC_1/PHOS_";
339 volpath += i;
340 volpath += "/PCPV_1";
341 if (gGeoManager->CheckPath(volpath.Data()))
342 isCPVpresent=1 ;
343 }
344
345 if(isCPVpresent){
346 nCPV = nEMC + geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() * nmod ;
347 }
348 else{
349 nCPV = nEMC ;
350 }
351
352 digits->Expand(nCPV) ;
353
354 //take all the inputs to add together and load the SDigits
355 TObjArray * sdigArray = new TObjArray(fInput) ;
356 sdigArray->AddAt(phosLoader->SDigits(), 0) ;
357
358 for(Int_t i = 1 ; i < fInput ; i++){
359 TString tempo(fEventNames[i]) ;
360 tempo += i ;
361 AliRunLoader* rl2 = AliRunLoader::GetRunLoader(tempo) ;
362 if(!rl2){
363 rl2 = AliRunLoader::Open(fInputFileNames[i], tempo) ;
364 if (!rl2) {
365 Fatal("AliPHOSDigitizer", "Could not find the Run Loader for %s - %s",fInputFileNames[i].Data(), tempo.Data()) ;
366 return ;
367 }
368 rl2->LoadHeader();
369 }
370 AliPHOSLoader * phosLoader2 = static_cast<AliPHOSLoader*>(rl2->GetLoader("PHOSLoader"));
371
372 if(fManager){
373 readEvent = static_cast<AliStream*>(fManager->GetInputStream(i))->GetCurrentEventNumber() ;
374 }
375 TClonesArray * digs ;
376 if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
377 AliInfo(Form("Adding event %d from input stream %d %s %s",
378 readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
379 rl2->GetEvent(readEvent) ;
380 phosLoader2->CleanDigits() ;
381 phosLoader2->LoadDigits("READ") ;
382 digs = phosLoader2->Digits() ;
383 toMakeNoise=0 ; //Do not add noise, it is already in stream
384 }
385 else{
386 AliInfo(Form("Adding event %d (SDigits) from input stream %d %s %s",
387 readEvent, i, fInputFileNames[i].Data(), tempo.Data())) ;
388 rl2->GetEvent(readEvent) ;
389 phosLoader2->CleanSDigits() ;
390 phosLoader2->LoadSDigits("READ") ;
391 digs = phosLoader2->SDigits() ;
392 }
393 sdigArray->AddAt(digs, i) ;
394 }
395
396 //Find the first crystal with signal
397 Int_t nextSig = 200000 ;
398 TClonesArray * sdigits ;
399 for(Int_t i = 0 ; i < fInput ; i++){
400 sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
401 if ( !sdigits->GetEntriesFast() )
402 continue ;
403 Int_t curNext = static_cast<AliPHOSDigit *>(sdigits->At(0))->GetId() ;
404 if(curNext < nextSig)
405 nextSig = curNext ;
406 }
407
408 TArrayI index(fInput) ;
409 index.Reset() ; //Set all indexes to zero
410
411 AliPHOSDigit * digit ;
412 AliPHOSDigit * curSDigit ;
413
414// TClonesArray * ticks = new TClonesArray("AliPHOSTick",1000) ;
415
416 //Put Noise contribution
417 Double_t apdNoise = 0. ;
418 if(toMakeNoise)
419 apdNoise = AliPHOSSimParam::GetInstance()->GetAPDNoise() ;
420
421 Int_t emcpermod=geom->GetNPhi()*geom->GetNZ();
422 Int_t idigit= 0;
423 for(Int_t imod=0; imod<5; imod++){
424 if(!isPresent[imod])
425 continue ;
426 Int_t firstAbsId=imod*emcpermod+1 ;
427 Int_t lastAbsId =(imod+1)*emcpermod ;
428 for(absID = firstAbsId ; absID <= lastAbsId ; absID++){
429 Float_t noise = gRandom->Gaus(0.,apdNoise) ;
430 new((*digits)[idigit]) AliPHOSDigit( -1, absID, noise, TimeOfNoise() ) ;
431 //look if we have to add signal?
432 digit = static_cast<AliPHOSDigit *>(digits->At(idigit)) ;
433 idigit++ ;
434
435 if(absID==nextSig){
436 //Add SDigits from all inputs
437// ticks->Clear() ;
438// Int_t contrib = 0 ;
439
440//New Timing model is necessary
441// Float_t a = digit->GetEnergy() ;
442// Float_t b = TMath::Abs( a / fTimeSignalLength) ;
443// //Mark the beginning of the signal
444// new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime(),0, b);
445// //Mark the end of the signal
446// new((*ticks)[contrib++]) AliPHOSTick(digit->GetTime()+fTimeSignalLength, -a, -b);
447
448// Calculate time as time of the largest digit
449 Float_t time = digit->GetTime() ;
450 Float_t eTime= digit->GetEnergy() ;
451
452 //loop over inputs
453 for(Int_t i = 0 ; i < fInput ; i++){
454 if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] ){
455 curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
456 if(AliPHOSSimParam::GetInstance()->IsStreamDigits(i)){ //This is Digits Stream
457 curSDigit->SetEnergy(Calibrate(curSDigit->GetEnergy(),curSDigit->GetId())) ;
458 curSDigit->SetTime(CalibrateT(curSDigit->GetTime(),curSDigit->GetId())) ;
459 }
460 }
461 else
462 curSDigit = 0 ;
463 //May be several digits will contribute from the same input
464 while(curSDigit && curSDigit->GetId() == absID){
465 //Shift primary to separate primaries belonging different inputs
466 Int_t primaryoffset ;
467 if(fManager)
468 primaryoffset = fManager->GetMask(i) ;
469 else
470 primaryoffset = 10000000*i ;
471 curSDigit->ShiftPrimary(primaryoffset) ;
472
473//New Timing model is necessary
474// a = curSDigit->GetEnergy() ;
475// b = a /fTimeSignalLength ;
476// new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime(),0, b);
477// new((*ticks)[contrib++]) AliPHOSTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b);
478 if(curSDigit->GetEnergy()>eTime){
479 eTime=curSDigit->GetEnergy() ;
480 time=curSDigit->GetTime() ;
481 }
482 *digit += *curSDigit ; //add energies
483
484 index[i]++ ;
485 if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
486 curSDigit = static_cast<AliPHOSDigit*>(static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i])) ;
487 else
488 curSDigit = 0 ;
489 }
490 }
491
492 //calculate and set time
493//New Timing model is necessary
494// Float_t time = FrontEdgeTime(ticks) ;
495 digit->SetTime(time) ;
496
497 //Find next signal module
498 nextSig = 200000 ;
499 for(Int_t i = 0 ; i < fInput ; i++){
500 sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
501 Int_t curNext = nextSig ;
502 if(sdigits->GetEntriesFast() > index[i] ){
503 curNext = static_cast<AliPHOSDigit *>(sdigits->At(index[i]))->GetId() ;
504 }
505 if(curNext < nextSig) nextSig = curNext ;
506 }
507 }
508 }
509 }
510
511 //distretize energy if necessary
512 if(AliPHOSSimParam::GetInstance()->IsEDigitizationOn()){
513 Float_t adcW=AliPHOSSimParam::GetInstance()->GetADCchannelW() ;
514 for(Int_t i = 0 ; i < nEMC ; i++){
515 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
516 digit->SetEnergy(adcW*ceil(digit->GetEnergy()/adcW)) ;
517 }
518 }
519 //Apply decalibration if necessary
520 for(Int_t i = 0 ; i < nEMC ; i++){
521 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
522 Decalibrate(digit) ;
523 }
524
525
526// ticks->Delete() ;
527// delete ticks ;
528
529 //Now CPV digits (different noise and no timing)
530 Int_t cpvpermod = geom->GetNumberOfCPVPadsZ() * geom->GetNumberOfCPVPadsPhi() ;
531 Int_t nEMCtotal=emcpermod*5 ;
532 Float_t cpvNoise = AliPHOSSimParam::GetInstance()->GetCPVNoise() ;
533 if(isCPVpresent){ //CPV is present in current geometry
534 for(Int_t imod=0; imod<5; imod++){ //module is present in current geometry
535 if(!isPresent[imod])
536 continue ;
537 Int_t firstAbsId=nEMCtotal+imod*cpvpermod+1 ;
538 Int_t lastAbsId =nEMCtotal+(imod+1)*cpvpermod ;
539 for(absID = firstAbsId; absID <= lastAbsId; absID++){
540 Float_t noise = gRandom->Gaus(0., cpvNoise) ;
541 new((*digits)[idigit]) AliPHOSDigit( -1,absID,noise, TimeOfNoise() ) ;
542 idigit++ ;
543 //look if we have to add signal?
544 if(absID==nextSig){
545 digit = static_cast<AliPHOSDigit *>(digits->At(idigit-1)) ;
546 //Add SDigits from all inputs
547 for(Int_t i = 0 ; i < fInput ; 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 //May be several digits will contribute from the same input
554 while(curSDigit && curSDigit->GetId() == absID){
555 //Shift primary to separate primaries belonging different inputs
556 Int_t primaryoffset ;
557 if(fManager)
558 primaryoffset = fManager->GetMask(i) ;
559 else
560 primaryoffset = 10000000*i ;
561 curSDigit->ShiftPrimary(primaryoffset) ;
562
563 //add energies
564 *digit += *curSDigit ;
565 index[i]++ ;
566 if( static_cast<TClonesArray *>(sdigArray->At(i))->GetEntriesFast() > index[i] )
567 curSDigit = static_cast<AliPHOSDigit*>( static_cast<TClonesArray *>(sdigArray->At(i))->At(index[i]) ) ;
568 else
569 curSDigit = 0 ;
570 }
571 }
572
573 //Find next signal module
574 nextSig = 200000 ;
575 for(Int_t i = 0 ; i < fInput ; i++){
576 sdigits = static_cast<TClonesArray *>(sdigArray->At(i)) ;
577 Int_t curNext = nextSig ;
578 if(sdigits->GetEntriesFast() > index[i] )
579 curNext = static_cast<AliPHOSDigit *>( sdigits->At(index[i]) )->GetId() ;
580 if(curNext < nextSig) nextSig = curNext ;
581 }
582
583 }
584 }
585 }
586 }
587
588 delete sdigArray ; //We should not delete its contents
589
590 Int_t relId[4];
591
592 //set amplitudes in bad channels to zero
593
594 for(Int_t i = 0 ; i <digits->GetEntriesFast(); i++){
595 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
596 geom->AbsToRelNumbering(digit->GetId(),relId);
597 if(relId[1] == 0) // Emc
598 if(fcdb->IsBadChannelEmc(relId[0],relId[3],relId[2])) digit->SetEnergy(0.);
599 }
600
601 //remove digits below thresholds
602 Float_t emcThreshold = AliPHOSSimParam::GetInstance()->GetEmcDigitsThreshold() ;
603 for(Int_t i = 0 ; i < nEMC ; i++){
604 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
605
606 if(digit->GetEnergy() < emcThreshold){
607 digits->RemoveAt(i) ;
608 continue ;
609 }
610
611 geom->AbsToRelNumbering(digit->GetId(),relId);
612
613 digit->SetEnergy(TMath::Ceil(digit->GetEnergy())-0.9999) ;
614
615 Float_t tres = TimeResolution(digit->GetEnergy()) ;
616 digit->SetTime(gRandom->Gaus(digit->GetTime(), tres) ) ;
617
618 fPulse->Reset();
619 fPulse->SetAmplitude(digit->GetEnergy()/
620 fcdb->GetADCchannelEmc(relId[0],relId[3],relId[2]));
621 fPulse->SetTZero(digit->GetTimeR());
622 fPulse->MakeSamples();
623 fPulse->GetSamples(fADCValuesHG, fADCValuesLG) ;
624 Int_t nSamples = fPulse->GetRawFormatTimeBins();
625 digit->SetALTROSamplesHG(nSamples,fADCValuesHG);
626 digit->SetALTROSamplesLG(nSamples,fADCValuesLG);
627 }
628
629 Float_t cpvDigitThreshold = AliPHOSSimParam::GetInstance()->GetCpvDigitsThreshold() ;
630 for(Int_t i = nEMC; i < nCPV ; i++){
631 if( static_cast<AliPHOSDigit*>(digits->At(i))->GetEnergy() < cpvDigitThreshold )
632 digits->RemoveAt(i) ;
633 }
634
635 digits->Compress() ;
636 Int_t ndigits = digits->GetEntriesFast() ;
637
638 //Set indexes in list of digits and make true digitization of the energy
639 for (Int_t i = 0 ; i < ndigits ; i++) {
640 digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
641 digit->SetIndexInList(i) ;
642 if(digit->GetId() > fEmcCrystals){ //digitize CPV only
643 digit->SetAmp(DigitizeCPV(digit->GetEnergy(),digit->GetId()) ) ;
644 }
645 }
646
647}
648//____________________________________________________________________________
649Float_t AliPHOSDigitizer::Calibrate(Float_t amp,Int_t absId){
650 //Apply calibration
651 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
652
653 //Determine rel.position of the cell absolute ID
654 Int_t relId[4];
655 geom->AbsToRelNumbering(absId,relId);
656 Int_t module=relId[0];
657 Int_t row =relId[2];
658 Int_t column=relId[3];
659 if(relId[1]==0){ //This Is EMC
660 Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
661 return amp*calibration ;
662 }
663 return 0 ;
664}
665//____________________________________________________________________________
666void AliPHOSDigitizer::Decalibrate(AliPHOSDigit *digit)
667{
668 // Decalibrate EMC digit, i.e. change its energy by a factor read from CDB
669
670 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
671
672 //Determine rel.position of the cell absolute ID
673 Int_t relId[4];
674 geom->AbsToRelNumbering(digit->GetId(),relId);
675 Int_t module=relId[0];
676 Int_t row =relId[2];
677 Int_t column=relId[3];
678 if(relId[1]==0){ //This Is EMC
679 Float_t calibration = fcdb->GetADCchannelEmc(module,column,row);
680 Float_t energy = digit->GetEnergy()/calibration;
681 digit->SetEnergy(energy); //Now digit measures E in ADC counts
682 Float_t time = digit->GetTime() ;
683 time-=fcdb->GetTimeShiftEmc(module,column,row);
684 digit->SetTime(time) ;
685 }
686}
687//____________________________________________________________________________
688Float_t AliPHOSDigitizer::CalibrateT(Float_t time,Int_t absId){
689 //Apply time calibration
690 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
691
692 //Determine rel.position of the cell absolute ID
693 Int_t relId[4];
694 geom->AbsToRelNumbering(absId,relId);
695 Int_t module=relId[0];
696 Int_t row =relId[2];
697 Int_t column=relId[3];
698 time += fcdb->GetTimeShiftEmc(module,column,row);
699 return time ;
700}
701//____________________________________________________________________________
702Int_t AliPHOSDigitizer::DigitizeCPV(Float_t charge, Int_t absId)
703{
704 // Returns digitized value of the CPV charge in a pad absId
705
706 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
707
708 //Determine rel.position of the cell absId
709 Int_t relId[4];
710 geom->AbsToRelNumbering(absId,relId);
711 Int_t module=relId[0];
712 Int_t row =relId[2];
713 Int_t column=relId[3];
714
715 Int_t channel = 0;
716
717 if(absId > fEmcCrystals){ //digitize CPV only
718
719 //reading calibration data for cell absId.
720 Float_t adcPedestalCpv = fcdb->GetADCpedestalCpv(module,column,row);
721 Float_t adcChanelCpv = fcdb->GetADCchannelCpv( module,column,row);
722
723 channel = (Int_t) TMath::Ceil((charge - adcPedestalCpv)/adcChanelCpv) ;
724 Int_t nMax = AliPHOSSimParam::GetInstance()->GetNADCcpv() ;
725 if(channel > nMax ) channel = nMax ;
726 }
727 return channel ;
728}
729
730//____________________________________________________________________________
731void AliPHOSDigitizer::Exec(Option_t *option)
732{
733 // Steering method to process digitization for events
734 // in the range from fFirstEvent to fLastEvent.
735 // This range is optionally set by SetEventRange().
736 // if fLastEvent=-1, then process events until the end.
737 // by default fLastEvent = fFirstEvent (process only one event)
738
739 if (!fInit) { // to prevent overwrite existing file
740 AliError(Form("Give a version name different from %s",
741 fEventFolderName.Data() )) ;
742 return ;
743 }
744
745 if (strstr(option,"print")) {
746 Print();
747 return ;
748 }
749
750 if(strstr(option,"tim"))
751 gBenchmark->Start("PHOSDigitizer");
752
753 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
754 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
755
756 // Post Digitizer to the white board
757 phosLoader->PostDigitizer(this) ;
758
759 if (fLastEvent == -1)
760 fLastEvent = rl->GetNumberOfEvents() - 1 ;
761 else if (fManager)
762 fLastEvent = fFirstEvent ;
763
764 Int_t nEvents = fLastEvent - fFirstEvent + 1;
765
766 Int_t ievent ;
767
768 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
769 fEventCounter++ ;
770
771 Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
772
773 WriteDigits() ;
774
775 if(strstr(option,"deb"))
776 PrintDigits(option);
777
778 //increment the total number of Digits per run
779 fDigitsInRun += phosLoader->Digits()->GetEntriesFast() ;
780 }
781
782 phosLoader->CleanDigitizer();
783
784 if(strstr(option,"tim")){
785 gBenchmark->Stop("PHOSDigitizer");
786 TString message ;
787 message = " took %f seconds for Digitizing %f seconds per event\n" ;
788 AliInfo(Form( message.Data(),
789 gBenchmark->GetCpuTime("PHOSDigitizer"),
790 gBenchmark->GetCpuTime("PHOSDigitizer")/nEvents ));
791 }
792}
793//____________________________________________________________________________
794Float_t AliPHOSDigitizer::TimeResolution(Float_t e){
795 //calculate TOF resolution using beam-test resutls
796 Float_t a=AliPHOSSimParam::GetInstance()->GetTOFa() ;
797 Float_t b=AliPHOSSimParam::GetInstance()->GetTOFb() ;
798 return TMath::Sqrt(a*a+b*b/e) ;
799}
800
801////____________________________________________________________________________
802//Float_t AliPHOSDigitizer::FrontEdgeTime(TClonesArray * ticks) const
803//{
804// // Returns the shortest time among all time ticks
805//
806// ticks->Sort() ; //Sort in accordance with times of ticks
807// TIter it(ticks) ;
808// AliPHOSTick * ctick = (AliPHOSTick *) it.Next() ;
809// Float_t time = ctick->CrossingTime(fTimeThreshold) ;
810//
811// AliPHOSTick * t ;
812// while((t=(AliPHOSTick*) it.Next())){
813// if(t->GetTime() < time) //This tick starts before crossing
814// *ctick+=*t ;
815// else
816// return time ;
817//
818// time = ctick->CrossingTime(fTimeThreshold) ;
819// }
820// return time ;
821//}
822
823//____________________________________________________________________________
824Bool_t AliPHOSDigitizer::Init()
825{
826 // Makes all memory allocations
827 fInit = kTRUE ;
828
829 AliPHOSGeometry *geom;
830 if (!(geom = AliPHOSGeometry::GetInstance()))
831 geom = AliPHOSGeometry::GetInstance("IHEP","");
832// const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
833
834 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
835
836 fFirstEvent = 0 ;
837 fLastEvent = fFirstEvent ;
838 if (fManager)
839 fInput = fManager->GetNinputs() ;
840 else
841 fInput = 1 ;
842
843 fInputFileNames = new TString[fInput] ;
844 fEventNames = new TString[fInput] ;
845 fInputFileNames[0] = GetTitle() ;
846 fEventNames[0] = fEventFolderName.Data() ;
847 Int_t index ;
848 for (index = 1 ; index < fInput ; index++) {
849 fInputFileNames[index] = static_cast<AliStream*>(fManager->GetInputStream(index))->GetFileName(0);
850 TString tempo = fManager->GetInputFolderName(index) ;
851 fEventNames[index] = tempo.Remove(tempo.Length()-1) ; // strip of the stream number added by fManager
852 }
853
854 //to prevent cleaning of this object while GetEvent is called
855 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
856 if(!rl){
857 rl = AliRunLoader::Open(GetTitle(), fEventFolderName) ;
858 }
859 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
860 phosLoader->GetDigitsDataLoader()->GetBaseTaskLoader()->SetDoNotReload(kTRUE);
861
862 return fInit ;
863}
864
865//____________________________________________________________________________
866void AliPHOSDigitizer::InitParameters()
867{
868 // Set initial parameters Digitizer
869
870 fDigitsInRun = 0 ;
871 SetEventRange(0,-1) ;
872 fPulse = new AliPHOSPulseGenerator();
873 fADCValuesLG = new Int_t[fPulse->GetRawFormatTimeBins()];
874 fADCValuesHG = new Int_t[fPulse->GetRawFormatTimeBins()];
875
876}
877
878//__________________________________________________________________
879void AliPHOSDigitizer::Print(const Option_t *)const
880{
881 // Print Digitizer's parameters
882 AliInfo(Form("\n------------------- %s -------------", GetName() )) ;
883 if( strcmp(fEventFolderName.Data(), "") != 0 ){
884 printf(" Writing Digits to branch with title %s\n", fEventFolderName.Data()) ;
885
886 Int_t nStreams ;
887 if (fManager)
888 nStreams = GetNInputStreams() ;
889 else
890 nStreams = fInput ;
891
892 Int_t index = 0 ;
893 for (index = 0 ; index < nStreams ; index++) {
894 TString tempo(fEventNames[index]) ;
895 tempo += index ;
896 printf ("Adding SDigits from %s \n", fInputFileNames[index].Data()) ;
897 }
898
899 // printf("\nWith following parameters:\n") ;
900 // printf(" Electronics noise in EMC (fPinNoise) = %f GeV\n", fPinNoise ) ;
901 // printf(" Threshold in EMC (fEMCDigitThreshold) = %f GeV\n", fEMCDigitThreshold ) ;
902 // printf(" Noise in CPV (fCPVNoise) = %f aux units\n", fCPVNoise ) ;
903 // printf(" Threshold in CPV (fCPVDigitThreshold) = %f aux units\n",fCPVDigitThreshold ) ;
904 printf(" ---------------------------------------------------\n") ;
905 }
906 else
907 AliInfo(Form("AliPHOSDigitizer not initialized" )) ;
908
909}
910
911//__________________________________________________________________
912 void AliPHOSDigitizer::PrintDigits(Option_t * option)
913{
914 // Print a table of digits
915
916
917 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
918 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
919 TClonesArray * digits = phosLoader->Digits() ;
920 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
921
922 AliInfo(Form("%d", digits->GetEntriesFast())) ;
923 printf("\nevent %d", gAlice->GetEvNumber()) ;
924 printf("\n Number of entries in Digits list %d", digits->GetEntriesFast() ) ;
925
926
927 if(strstr(option,"all")||strstr(option,"EMC")){
928 //loop over digits
929 AliPHOSDigit * digit;
930 printf("\nEMC 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()) &&
935 (static_cast<AliPHOSDigit *>(digits->At(index))->GetId() <= maxEmc) ; index++) {
936 digit = (AliPHOSDigit * ) digits->At(index) ;
937 if(digit->GetNprimary() == 0)
938 continue;
939// printf("%6d %8d %6.5e %4d %2d :",
940// digit->GetId(), digit->GetAmp(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ; // YVK
941 printf("%6d %.4f %6.5e %4d %2d :",
942 digit->GetId(), digit->GetEnergy(), digit->GetTime(), digit->GetIndexInList(), digit->GetNprimary()) ;
943 Int_t iprimary;
944 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
945 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
946 }
947 printf("\n") ;
948 }
949 }
950
951 if(strstr(option,"all")||strstr(option,"CPV")){
952
953 //loop over CPV digits
954 AliPHOSDigit * digit;
955 printf("\nCPV digits (with primaries):\n") ;
956 printf("\n Id Amplitude Time Index Nprim: Primaries list \n") ;
957 Int_t maxEmc = geom->GetNModules()*geom->GetNCristalsInModule() ;
958 Int_t index ;
959 for (index = 0 ; index < digits->GetEntriesFast(); index++) {
960 digit = (AliPHOSDigit * ) digits->At(index) ;
961 if(digit->GetId() > maxEmc){
962 printf("%6d %8d %4d %2d :",
963 digit->GetId(), digit->GetAmp(), digit->GetIndexInList(), digit->GetNprimary()) ;
964 Int_t iprimary;
965 for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++) {
966 printf("%d ",digit->GetPrimary(iprimary+1) ) ;
967 }
968 printf("\n") ;
969 }
970 }
971 }
972
973}
974
975//__________________________________________________________________
976Float_t AliPHOSDigitizer::TimeOfNoise(void) const
977{ // Calculates the time signal generated by noise
978 //PH Info("TimeOfNoise", "Change me") ;
979 return gRandom->Rndm() * 1.28E-5;
980}
981
982//__________________________________________________________________
983void AliPHOSDigitizer::Unload()
984{
985
986 Int_t i ;
987 for(i = 1 ; i < fInput ; i++){
988 TString tempo(fEventNames[i]) ;
989 tempo += i ;
990 AliRunLoader* rl = AliRunLoader::GetRunLoader(tempo) ;
991 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
992 phosLoader->UnloadSDigits() ;
993 }
994
995 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
996 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
997 phosLoader->UnloadDigits() ;
998}
999
1000//____________________________________________________________________________
1001void AliPHOSDigitizer::WriteDigits()
1002{
1003
1004 // Makes TreeD in the output file.
1005 // Check if branch already exists:
1006 // if yes, exit without writing: ROOT TTree does not support overwriting/updating of
1007 // already existing branches.
1008 // else creates branch with Digits, named "PHOS", title "...",
1009 // and branch "AliPHOSDigitizer", with the same title to keep all the parameters
1010 // and names of files, from which digits are made.
1011
1012 AliRunLoader* rl = AliRunLoader::GetRunLoader(fEventFolderName) ;
1013 AliPHOSLoader * phosLoader = static_cast<AliPHOSLoader*>(rl->GetLoader("PHOSLoader"));
1014
1015 const TClonesArray * digits = phosLoader->Digits() ;
1016 TTree * treeD = phosLoader->TreeD();
1017 if(!treeD){
1018 phosLoader->MakeTree("D");
1019 treeD = phosLoader->TreeD();
1020 }
1021
1022 // -- create Digits branch
1023 Int_t bufferSize = 32000 ;
1024 TBranch * digitsBranch = treeD->Branch("PHOS","TClonesArray",&digits,bufferSize);
1025 digitsBranch->SetTitle(fEventFolderName);
1026 digitsBranch->Fill() ;
1027
1028 phosLoader->WriteDigits("OVERWRITE");
1029 phosLoader->WriteDigitizer("OVERWRITE");
1030
1031 Unload() ;
1032
1033}