Default arguments set only in the header file
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALDigitizer.cxx
1 /**************************************************************************
2
3
4  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
5
6
7  *                                                                        *
8
9
10  * Author: The ALICE Off-line Project.                                    *
11
12
13  * Contributors are mentioned in the code where appropriate.              *
14
15
16  *                                                                        *
17
18
19  * Permission to use, copy, modify and distribute this software and its   *
20
21
22  * documentation strictly for non-commercial purposes is hereby granted   *
23
24
25  * without fee, provided that the above copyright notice appears in all   *
26
27
28  * copies and that both the copyright notice and this permission notice   *
29
30
31  * appear in the supporting documentation. The authors make no claims     *
32
33
34  * about the suitability of this software for any purpose. It is          *
35
36
37  * provided "as is" without express or implied warranty.                  *
38
39
40  **************************************************************************/
41
42
43
44
45
46 /* $Id$ */
47
48
49
50
51
52 //_________________________________________________________________________
53
54
55 // 
56
57
58 //////////////////////////////////////////////////////////////////////////////
59
60
61 // Class performs digitization of Summable digits 
62
63
64 //  
65
66
67 // In addition it performs mixing of summable digits from different events.
68
69
70 //
71
72
73 // For each event two branches are created in TreeD:
74
75
76 //   "EMCAL" - list of digits
77
78
79 //   "AliEMCALDigitizer" - AliEMCALDigitizer with all parameters used in digitization
80
81
82 //
83
84
85 // Note, that one cset title for new digits branch, and repeat digitization with
86
87
88 // another set of parameters.
89
90
91 //
92
93
94 // Examples of use:
95
96
97 // root[0] AliEMCALDigitizer * d = new AliEMCALDigitizer() ;
98
99
100 // root[1] d->ExecuteTask()             
101
102
103 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
104
105
106 //                       //Digitizes SDigitis in all events found in file galice.root 
107
108
109 //
110
111
112 // root[2] AliEMCALDigitizer * d1 = new AliEMCALDigitizer("galice1.root") ;  
113
114
115 //                       // Will read sdigits from galice1.root
116
117
118 // root[3] d1->MixWith("galice2.root")       
119
120
121 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
122
123
124 //                       // Reads another portion of sdigits from galice2.root
125
126
127 // root[3] d1->MixWith("galice3.root")       
128
129
130 //                       // Reads another portion of sdigits from galice3.root
131
132
133 // root[4] d->ExecuteTask("deb timing")    
134
135
136 //                       // Reads SDigits from files galice1.root, galice2.root ....
137
138
139 //                       // mixes them and stores produced Digits in file galice1.root          
140
141
142 //                       // deb - prints number of produced digits
143
144
145 //                       // deb all - prints list of produced digits
146
147
148 //                       // timing  - prints time used for digitization
149
150
151 ////////////////////////////////////////////////////////////////////////////////////
152
153
154 //
155
156
157 //*-- Author: Sahal Yacoob (LBL)
158
159
160 // based on : AliEMCALDigitizer
161
162
163 //_________________________________________________________________________________
164
165
166
167
168
169 // --- ROOT system ---
170
171
172 #include "TFile.h"
173
174
175 #include "TTree.h"
176
177
178 #include "TSystem.h"
179
180
181 #include "TROOT.h"
182
183
184 #include "TFolder.h"
185
186
187 #include "TObjString.h"
188
189
190 #include "TBenchmark.h"
191
192
193 // --- Standard library ---
194
195
196 #include <iomanip.h>
197
198
199
200
201
202 // --- AliRoot header files ---
203
204
205
206
207
208 #include "AliRun.h"
209
210
211 #include "AliEMCALDigit.h"
212
213
214 #include "AliEMCALHit.h"
215
216
217 #include "AliEMCALTick.h"
218
219
220 #include "AliEMCALv1.h"
221
222
223 #include "AliEMCALDigitizer.h"
224
225
226 #include "AliEMCALSDigitizer.h"
227
228
229 #include "AliEMCALGeometry.h"
230
231
232 #include "AliEMCALGetter.h"
233
234
235 #include "AliRunDigitizer.h"
236
237
238 ClassImp(AliEMCALDigitizer)
239
240
241
242
243
244
245
246
247 //____________________________________________________________________________ 
248
249
250   AliEMCALDigitizer::AliEMCALDigitizer()
251
252
253 {
254
255
256   // ctor
257
258
259
260
261
262   fSDigitizer = 0 ;
263
264
265   fNinputs = 0 ;
266
267
268   fPinNoise = 0.0 ;
269
270
271   fTowerDigitThreshold = 0.0 ;
272
273
274   fTimeResolution     = 0. ;
275
276
277   fTimeSignalLength   = 0. ;
278
279
280   fPreShowerDigitThreshold = 0. ;
281
282
283   fADCchannelTower = 0.0;      // width of one ADC channel in GeV
284
285
286   fADCpedestalTower = 0. ;      // pedestal of ADC
287
288
289   fNADCTower = 0;  // number of channels in Tower ADC
290
291
292
293
294
295   fADCchannelPreSho  = 0.0;          // width of one ADC channel in Pre Shower
296
297
298   fADCpedestalPreSho = 0.0 ;         // pedestal of ADC
299
300
301   fNADCPreSho = 0;      // number of channels in Pre Shower ADC
302
303
304   fTimeThreshold = 0.0; //Means 1 MeV in terms of SDigits amplitude
305
306
307   fManager = 0 ;
308
309
310
311
312
313
314
315
316
317
318
319 }
320
321
322 //____________________________________________________________________________ 
323
324
325 Bool_t AliEMCALDigitizer::Init()
326
327
328 {
329
330
331   // Makes all memory allocations
332
333
334
335
336
337   fSDigitizer = 0 ;
338
339
340   fNinputs = 1 ;
341
342
343   fPinNoise = 0.00001 ;
344
345
346   fTowerDigitThreshold = 0.001 ;
347
348
349   fTimeResolution     = 0.5e-9 ;
350
351
352   fTimeSignalLength   = 1.0e-9 ;
353
354
355   fPreShowerDigitThreshold = fTowerDigitThreshold/25. ;
356
357
358   fInitialized = kFALSE ;
359
360
361   fADCchannelTower = 0.000220;       // width of one ADC channel in GeV
362
363
364   fADCpedestalTower = 0.005 ;      // GeV
365
366
367   fNADCTower = (Int_t) TMath::Power(2,16) ;  // number of channels in Tower ADC
368
369
370
371
372
373   fADCchannelPreSho = 0.0000300;          // width of one ADC channel in Pre Shower
374
375
376   fADCpedestalPreSho = 0.005 ;         // 
377
378
379   fNADCPreSho = (Int_t) TMath::Power(2,12);      // number of channels in Pre ShowerADC
380
381
382
383
384
385   fTimeThreshold = 0.001*10000000 ; //Means 1 MeV in terms of SDigits amplitude
386
387
388  
389
390
391
392
393
394
395
396
397 if(fManager)
398
399
400     SetTitle("aliroot") ;
401
402
403   else if (strcmp(GetTitle(),"")==0) 
404
405
406    SetTitle("galice.root") ;
407
408
409
410
411
412   if( strcmp(GetName(), "") == 0 )
413
414
415     SetName("Default") ;
416
417
418   
419
420
421   AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), GetName(), "update") ; 
422
423
424   if ( gime == 0 ) {
425
426
427     cerr << "ERROR: AliEMCALDigitizer::Init -> Could not obtain the Getter object !" << endl ; 
428
429
430     return kFALSE;
431
432
433   } 
434
435
436   
437
438
439   //const AliEMCALGeometry * geom = gime->EMCALGeometry() ;
440
441
442   //fEmcCrystals = geom->GetNModules() *  geom->GetNCristalsInModule() ;
443
444
445   
446
447
448   // Post Digits to the white board
449
450
451   gime->PostDigits(GetName() ) ;   
452
453
454   
455
456
457   // Post Digitizer to the white board
458
459
460   gime->PostDigitizer(this) ;
461
462
463   
464
465
466   //Mark that we will use current header file
467
468
469   if(!fManager){
470
471
472     gime->PostSDigits(GetName(),GetTitle()) ;
473
474
475     gime->PostSDigitizer(GetName(),GetTitle()) ;
476
477
478   }
479
480
481   return kTRUE ;
482
483
484   
485
486
487  
488
489
490   
491
492
493 }
494
495
496
497
498
499 //____________________________________________________________________________ 
500
501
502 AliEMCALDigitizer::AliEMCALDigitizer(const char *headerFile,const char *name)
503
504
505 {
506
507
508    SetName(name) ;
509
510
511   SetTitle(headerFile) ;
512
513
514   fManager = 0 ;                     // We work in the standalong mode
515
516
517   Init() ;
518
519
520   
521
522
523
524
525
526   
527
528
529 }
530
531
532 //____________________________________________________________________________ 
533
534
535 AliEMCALDigitizer::AliEMCALDigitizer(AliRunDigitizer * ard):AliDigitizer(ard)
536
537
538 {
539
540
541   // ctor
542
543
544   SetName("");     //Will call init in the digitizing
545
546
547   SetTitle("aliroot") ;  
548
549
550 }
551
552
553
554
555
556 //____________________________________________________________________________ 
557
558
559   AliEMCALDigitizer::~AliEMCALDigitizer()
560
561
562 {
563
564
565   // dtor
566
567
568
569
570
571 }
572
573
574 //____________________________________________________________________________
575
576
577 void AliEMCALDigitizer::Reset() { 
578
579
580   //sets current event number to the first simulated event
581
582
583 if( strcmp(GetName(), "") == 0 )
584
585
586   Init() ;
587
588
589
590
591
592       // Int_t inputs ;
593
594
595       // for(inputs = 0; inputs < fNinputs ;inputs++)
596
597
598       //  fIevent->AddAt(-1, inputs ) ;
599
600
601   
602
603
604 }
605
606
607
608
609
610 //____________________________________________________________________________
611
612
613 void AliEMCALDigitizer::Digitize(const Int_t event) { 
614
615
616
617
618
619   // Makes the digitization of the collected summable digits
620
621
622   // for this it first creates the array of all EMCAL modules
623
624
625   // filled with noise (different for EMC, CPV and PPSD) and
626
627
628   // after that adds contributions from SDigits. This design 
629
630
631   // helps to avoid scanning over the list of digits to add 
632
633
634   // contribution of any new SDigit.
635
636
637
638
639
640   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
641
642
643   TClonesArray * digits = gime->Digits(GetName()) ; 
644
645
646   
647
648
649   digits->Clear() ;
650
651
652
653
654
655     const AliEMCALGeometry *geom = gime->EMCALGeometry() ; 
656
657
658
659
660
661
662
663
664   //Making digits with noise, first EMC
665
666
667   Int_t nEMC = 2*geom->GetNPhi()*geom->GetNZ();
668
669
670   Int_t absID ;
671
672
673   TString name      =  geom->GetName() ;
674
675
676
677
678
679  // get first the sdigitizer from the tasks list (must have same name as the digitizer)
680
681
682   const AliEMCALSDigitizer * sDigitizer = gime->SDigitizer(GetName()); 
683
684
685   if ( !sDigitizer) {
686
687
688     cerr << "ERROR: AliEMCALDigitizer::Digitize -> SDigitizer with name " << GetName() << " not found " << endl ; 
689
690
691     abort() ; 
692
693
694   }
695
696
697
698
699
700 // loop through the sdigits posted to the White Board and add them to the noise
701
702
703   TCollection * folderslist = gime->SDigitsFolder()->GetListOfFolders() ; 
704
705
706   TIter next(folderslist) ; 
707
708
709   TFolder * folder = 0 ; 
710
711
712   TClonesArray * sdigits = 0 ;
713
714
715   Int_t input = 0 ;
716
717
718   TObjArray * sdigArray = new TObjArray(2) ;
719
720
721   while ( (folder = (TFolder*)next()) ) 
722
723
724     if ( (sdigits = (TClonesArray*)folder->FindObject(GetName()) ) ) {
725
726
727       cout << "INFO: AliEMCALDigitizer::Digitize -> Adding SDigits " 
728
729
730            << GetName() << " from " << folder->GetName() << endl ; 
731
732
733       sdigArray->AddAt(sdigits, input) ;
734
735
736       input++ ;
737
738
739     }
740
741
742
743
744
745   //Find the first tower with signal
746
747
748   Int_t nextSig = 200000 ; 
749
750
751   Int_t i;
752
753
754   for(i=0; i<input; i++){
755
756
757     sdigits = (TClonesArray *)sdigArray->At(i) ;
758
759
760     if ( !sdigits->GetEntriesFast() )
761
762
763       continue ; 
764
765
766     Int_t curNext = ((AliEMCALDigit *)sdigits->At(0))->GetId() ;
767
768
769      if(curNext < nextSig) 
770
771
772        nextSig = curNext ;
773
774
775   }
776
777
778
779
780
781   TArrayI index(input) ;
782
783
784   index.Reset() ;  //Set all indexes to zero
785
786
787
788
789
790   AliEMCALDigit * digit = 0 ;
791
792
793   AliEMCALDigit * curSDigit = 0 ;
794
795
796
797
798
799   TClonesArray * ticks = new TClonesArray("AliEMCALTick",1000) ;
800
801
802
803
804
805   //Put Noise contribution
806
807
808   for(absID = 1; absID <= nEMC; absID++){
809
810
811     Float_t noise = gRandom->Gaus(0., fPinNoise); 
812
813
814     new((*digits)[absID-1]) AliEMCALDigit( -1, -1, absID,sDigitizer->Digitize(noise), TimeOfNoise() ) ;
815
816
817     //look if we have to add signal?
818
819
820     if(absID==nextSig){
821
822
823       //Add SDigits from all inputs 
824
825
826       digit = (AliEMCALDigit *) digits->At(absID-1) ;
827
828
829       
830
831
832       ticks->Clear() ;
833
834
835       Int_t contrib = 0 ;
836
837
838       Float_t a = digit->GetAmp() ;
839
840
841       Float_t b = TMath::Abs( a /fTimeSignalLength) ;
842
843
844       //Mark the beginnign of the signal
845
846
847       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime(),0, b);  
848
849
850       //Mark the end of the ignal     
851
852
853       new((*ticks)[contrib++]) AliEMCALTick(digit->GetTime()+fTimeSignalLength, -a, -b);
854
855
856       
857
858
859  // loop over input
860
861
862   
863
864
865       for(i = 0; i< input ; i++){  //loop over (possible) merge sources
866
867
868         if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
869
870
871           curSDigit = (AliEMCALDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;        
872
873
874         else
875
876
877           curSDigit = 0 ;
878
879
880         //May be several digits will contribute from the same input
881
882
883         while(curSDigit && curSDigit->GetId() == absID){           
884
885
886           //Shift primary to separate primaries belonging different inputs
887
888
889           Int_t primaryoffset ;
890
891
892           if(fManager)
893
894
895             primaryoffset = fManager->GetMask(i) ; 
896
897
898           else
899
900
901             primaryoffset = i ;
902
903
904           curSDigit->ShiftPrimary(primaryoffset) ;
905
906
907           
908
909
910           a = curSDigit->GetAmp() ;
911
912
913           b = a /fTimeSignalLength ;
914
915
916           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime(),0, b);  
917
918
919           new((*ticks)[contrib++]) AliEMCALTick(curSDigit->GetTime()+fTimeSignalLength, -a, -b); 
920
921
922
923
924
925           *digit = *digit + *curSDigit ;  //add energies
926
927
928
929
930
931           index[i]++ ;
932
933
934           if(((TClonesArray *)sdigArray->At(i))->GetEntriesFast() > index[i] )
935
936
937             curSDigit = (AliEMCALDigit*)((TClonesArray *)sdigArray->At(i))->At(index[i]) ;      
938
939
940           else
941
942
943             curSDigit = 0 ;
944
945
946         }
947
948
949       }
950
951
952
953
954
955 //calculate and set time
956
957
958       Float_t time = FrontEdgeTime(ticks) ;
959
960
961       digit->SetTime(time) ;
962
963
964
965
966
967       //Find next signal module
968
969
970       nextSig = 200000 ;
971
972
973       for(i=0; i<input; i++){
974
975
976         sdigits = ((TClonesArray *)sdigArray->At(i)) ;
977
978
979         Int_t curNext = nextSig ;
980
981
982         if(sdigits->GetEntriesFast() > index[i] ){
983
984
985           curNext = ((AliEMCALDigit *) sdigits->At(index[i]))->GetId() ;
986
987
988           
989
990
991         }
992
993
994         if(curNext < nextSig) nextSig = curNext ;
995
996
997       }
998
999
1000     }
1001
1002
1003   }
1004
1005
1006   
1007
1008
1009   ticks->Delete() ;
1010
1011
1012   delete ticks ;
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027   //remove digits below thresholds
1028
1029
1030   for(absID = 0; absID < nEMC/2 ; absID++){
1031
1032
1033     if(sDigitizer->Calibrate(((AliEMCALDigit*)digits->At(absID))->GetAmp()) < fTowerDigitThreshold)
1034
1035
1036       digits->RemoveAt(absID) ;
1037
1038
1039     else
1040
1041
1042       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
1043
1044
1045   }
1046
1047
1048   
1049
1050
1051   for(absID = nEMC/2; absID < nEMC ; absID++){
1052
1053
1054     if(sDigitizer->Calibrate(((AliEMCALDigit*)digits->At(absID))->GetAmp()) < fPreShowerDigitThreshold)
1055
1056
1057       digits->RemoveAt(absID) ;
1058
1059
1060     else
1061
1062
1063       digit->SetTime(gRandom->Gaus(digit->GetTime(),fTimeResolution) ) ;
1064
1065
1066   }
1067
1068
1069   
1070
1071
1072   digits->Compress() ;  
1073
1074
1075   
1076
1077
1078   Int_t ndigits = digits->GetEntriesFast() ;
1079
1080
1081   
1082
1083
1084   digits->Expand(ndigits) ;
1085
1086
1087   
1088
1089
1090   
1091
1092
1093   //Set indexes in list of digits
1094
1095
1096   //Int_t i ;
1097
1098
1099  for (i = 0 ; i < ndigits ; i++) { 
1100
1101
1102     AliEMCALDigit * digit = (AliEMCALDigit *) digits->At(i) ; 
1103
1104
1105     digit->SetIndexInList(i) ; 
1106
1107
1108     Float_t energy = sDigitizer->Calibrate(digit->GetAmp()) ;
1109
1110
1111     digit->SetAmp(DigitizeEnergy(energy,digit->GetId()) ) ;
1112
1113
1114   }
1115
1116
1117 }
1118
1119
1120
1121
1122
1123 //____________________________________________________________________________
1124
1125
1126
1127
1128
1129 Int_t AliEMCALDigitizer::DigitizeEnergy(Float_t energy, Int_t absId)
1130
1131
1132
1133
1134
1135   Int_t channel = -999;
1136
1137
1138   Int_t nphi = AliEMCALGetter::GetInstance()->EMCALGeometry()->GetNPhi() ; 
1139
1140
1141   Int_t nz   = AliEMCALGetter::GetInstance()->EMCALGeometry()->GetNZ() ;
1142
1143
1144   
1145
1146
1147   if(absId <= nphi*nz){  //digitize as tower
1148
1149
1150     channel = static_cast<Int_t> (TMath::Ceil( (energy + fADCpedestalTower)/fADCchannelTower ))  ;
1151
1152
1153   if(channel > fNADCTower ) 
1154
1155
1156     channel =  fNADCTower ;
1157
1158
1159   } else {
1160
1161
1162     channel =  static_cast<Int_t>(TMath::Ceil( (energy + fADCpedestalPreSho)/fADCchannelPreSho ))  ;
1163
1164
1165   if(channel > fNADCPreSho ) 
1166
1167
1168     channel =  fNADCPreSho ;
1169
1170
1171   }
1172
1173
1174   
1175
1176
1177   return channel ;
1178
1179
1180 }
1181
1182
1183
1184
1185
1186 //____________________________________________________________________________
1187
1188
1189 void AliEMCALDigitizer::Exec(Option_t *option) { 
1190
1191
1192   // Managing method
1193
1194
1195 if(strcmp(GetName(), "") == 0 )   
1196
1197
1198     Init() ;
1199
1200
1201   
1202
1203
1204   if (strstr(option,"print")) {
1205
1206
1207     Print("");
1208
1209
1210     return ; 
1211
1212
1213   }
1214
1215
1216   
1217
1218
1219   if(strstr(option,"tim"))
1220
1221
1222     gBenchmark->Start("EMCALDigitizer");
1223
1224
1225
1226
1227
1228   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
1229
1230
1231   
1232
1233
1234   Int_t nevents ;
1235
1236
1237   
1238
1239
1240   TTree * treeD ;
1241
1242
1243   
1244
1245
1246   if(fManager){
1247
1248
1249     treeD = fManager->GetTreeD() ;
1250
1251
1252     nevents = 1 ;    // Will process only one event
1253
1254
1255   }
1256
1257
1258   else {
1259
1260
1261     gAlice->GetEvent(0) ;
1262
1263
1264     nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
1265
1266
1267     treeD=gAlice->TreeD() ;
1268
1269
1270   }
1271
1272
1273   if(treeD == 0 ){
1274
1275
1276     cerr << " AliEMCALDigitizer :: Can not find TreeD " << endl ;
1277
1278
1279     return ;
1280
1281
1282   }
1283
1284
1285
1286
1287
1288   //Check, if this branch already exits
1289
1290
1291   TObjArray * lob = (TObjArray*)treeD->GetListOfBranches() ;
1292
1293
1294   TIter next(lob) ; 
1295
1296
1297   TBranch * branch = 0 ;  
1298
1299
1300   Bool_t emcalfound = kFALSE, digitizerfound = kFALSE ; 
1301
1302
1303   
1304
1305
1306   while ( (branch = (TBranch*)next()) && (!emcalfound || !digitizerfound) ) {
1307
1308
1309     if ( (strcmp(branch->GetName(), "EMCAL")==0) && 
1310
1311
1312          (strcmp(branch->GetTitle(), GetName())==0) ) 
1313
1314
1315       emcalfound = kTRUE ;
1316
1317
1318     
1319
1320
1321     else if ( (strcmp(branch->GetName(), "AliEMCALDigitizer")==0) && 
1322
1323
1324               (strcmp(branch->GetTitle(), GetName())==0) ) 
1325
1326
1327       digitizerfound = kTRUE ; 
1328
1329
1330   }
1331
1332
1333
1334
1335
1336   if ( emcalfound ) {
1337
1338
1339     cerr << "WARNING: AliEMCALDigitizer -> Digits branch with name " << GetName() 
1340
1341
1342          << " already exits" << endl ;
1343
1344
1345     return ; 
1346
1347
1348   }   
1349
1350
1351   if ( digitizerfound ) {
1352
1353
1354     cerr << "WARNING: AliEMCALDigitizer -> Digitizer branch with name " << GetName() 
1355
1356
1357          << " already exits" << endl ;
1358
1359
1360     return ; 
1361
1362
1363   }   
1364
1365
1366
1367
1368
1369   Int_t ievent ;
1370
1371
1372
1373
1374
1375   for(ievent = 0; ievent < nevents; ievent++){
1376
1377
1378     
1379
1380
1381     if(fManager){
1382
1383
1384       Int_t input ;
1385
1386
1387       for(input = 0 ; input < fManager->GetNinputs(); input ++){
1388
1389
1390         TTree * treeS = fManager->GetInputTreeS(input) ;
1391
1392
1393         if(!treeS){
1394
1395
1396           cerr << "AliEMCALDigitizer -> No Input " << endl ;
1397
1398
1399           return ;
1400
1401
1402         }
1403
1404
1405         gime->ReadTreeS(treeS,input) ;
1406
1407
1408       }
1409
1410
1411     }
1412
1413
1414     else
1415
1416
1417       gime->Event(ievent,"S") ;
1418
1419
1420     
1421
1422
1423     Digitize(ievent) ; //Add prepared SDigits to digits and add the noise
1424
1425
1426     
1427
1428
1429     WriteDigits(ievent) ;
1430
1431
1432     
1433
1434
1435     if(strstr(option,"deb"))
1436
1437
1438       PrintDigits(option);
1439
1440
1441     
1442
1443
1444     //increment the total number of Digits per run 
1445
1446
1447     fDigitsInRun += gime->Digits()->GetEntriesFast() ;  
1448
1449
1450   }
1451
1452
1453   
1454
1455
1456   if(strstr(option,"tim")){
1457
1458
1459     gBenchmark->Stop("EMCALDigitizer");
1460
1461
1462     cout << "AliEMCALDigitizer:" << endl ;
1463
1464
1465     cout << "  took " << gBenchmark->GetCpuTime("EMCALDigitizer") << " seconds for Digitizing " 
1466
1467
1468          <<  gBenchmark->GetCpuTime("EMCALDigitizer")/nevents << " seconds per event " << endl ;
1469
1470
1471     cout << endl ;
1472
1473
1474   }
1475
1476
1477   
1478
1479
1480 }
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492 //__________________________________________________________________
1493
1494
1495 void AliEMCALDigitizer::MixWith(char* headerFile){
1496
1497
1498   // Alows produce digits by superimposing background and signal event.
1499
1500
1501   // It is assumed, that headers file with SIGNAL events is opened in 
1502
1503
1504   // constructor, and now we set the BACKGROUND event, with which we 
1505
1506
1507   // will mix. Thus we avoid writing (changing) huge and expencive 
1508
1509
1510   // backgound files: all output will be writen into SIGNAL, i.e. 
1511
1512
1513   // opened in constructor file. 
1514
1515
1516   //
1517
1518
1519   // One can open as many files to mix with as one wants.
1520
1521
1522
1523
1524
1525 if( strcmp(GetName(), "") == 0 )
1526
1527
1528     Init() ;
1529
1530
1531   
1532
1533
1534  if(fManager){
1535
1536
1537     cout << "Can not use this method under AliRunDigitizer " << endl ;
1538
1539
1540     return ;
1541
1542
1543   } // check if the specified SDigits do not already exist on the White Board:
1544
1545
1546   // //Folders/RunMC/Event/Data/EMCAL/SDigits/headerFile/sdigitsname
1547
1548
1549
1550
1551
1552   TString path = "Folders/RunMC/Event/Data/EMCAL/SDigits" ; 
1553
1554
1555   path += headerFile ; 
1556
1557
1558   path += "/" ; 
1559
1560
1561   path += GetName() ;
1562
1563
1564   if ( gROOT->FindObjectAny(path.Data()) ) {
1565
1566
1567     cerr << "WARNING: AliEMCALDigitizer::MixWith -> Entry already exists, do not add" << endl ;
1568
1569
1570     return;
1571
1572
1573   }
1574
1575
1576
1577
1578
1579   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ;
1580
1581
1582   gime->PostSDigits(GetName(),headerFile) ;
1583
1584
1585   
1586
1587
1588   // check if the requested file is already open or exist and if SDigits Branch exist
1589
1590
1591   TFile * file = (TFile*)gROOT->FindObject(headerFile); 
1592
1593
1594   if ( !file ) { 
1595
1596
1597     file = new TFile(headerFile, "READ") ; 
1598
1599
1600     if (!file) { 
1601
1602
1603       cerr << "ERROR: AliEMCALDigitizer::MixWith -> File " << headerFile << " does not exist!" << endl ; 
1604
1605
1606       return ; 
1607
1608
1609     }
1610
1611
1612   }
1613
1614
1615   
1616
1617
1618 }
1619
1620
1621  
1622
1623
1624 //__________________________________________________________________
1625
1626
1627 void AliEMCALDigitizer::Print(Option_t* option)const {
1628
1629
1630   if( strcmp(GetName(), "") != 0) {
1631
1632
1633     
1634
1635
1636     cout << "------------------- "<< GetName() << " -------------" << endl ;
1637
1638
1639     cout << "Digitizing sDigits from file(s): " <<endl ;
1640
1641
1642     
1643
1644
1645   TCollection * folderslist = ((TFolder*)gROOT->FindObjectAny("Folders/RunMC/Event/Data/EMCAL/SDigits"))->GetListOfFolders() ; 
1646
1647
1648     TIter next(folderslist) ; 
1649
1650
1651     TFolder * folder = 0 ;
1652
1653
1654     while ( (folder = (TFolder*)next()) ) 
1655
1656
1657       if ( folder->FindObject(GetName())  ) 
1658
1659
1660         {
1661
1662
1663         cout << "Adding SDigits " << GetName() << " from " << folder->GetName() << endl ; 
1664
1665
1666       cout << endl ;
1667
1668
1669       cout << "Writing digits to " << GetTitle() << endl ;
1670
1671
1672       
1673
1674
1675       cout << endl ;
1676
1677
1678       cout << "With following parameters: " << endl ;
1679
1680
1681       cout << "     Electronics noise in EMC (fPinNoise) = " << fPinNoise << endl ;
1682
1683
1684       cout << "  Threshold  in EMC  (fTowerDigitThreshold) = " << fTowerDigitThreshold  << endl;
1685
1686
1687       cout << "  Threshold  in PreShower  (fPreShowerDigitThreshold) = " << fPreShowerDigitThreshold  << endl ; ;
1688
1689
1690       cout << "---------------------------------------------------" << endl ;
1691
1692
1693     }
1694
1695
1696     else
1697
1698
1699       cout << "AliEMCALDigitizer not initialized " << endl ;
1700
1701
1702     }
1703
1704
1705 }
1706
1707
1708
1709
1710
1711 //__________________________________________________________________
1712
1713
1714 void AliEMCALDigitizer::PrintDigits(Option_t * option){
1715
1716
1717     
1718
1719
1720   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
1721
1722
1723   TClonesArray * fDigits = gime->Digits() ;
1724
1725
1726
1727
1728
1729   cout << "AliEMCALDigitiser:"<< endl ;
1730
1731
1732   cout << "       Number of entries in Digits list " << fDigits->GetEntriesFast() << endl ;
1733
1734
1735   cout << endl ;
1736
1737
1738   if(strstr(option,"all")){
1739
1740
1741     
1742
1743
1744     //loop over digits
1745
1746
1747     AliEMCALDigit * digit;
1748
1749
1750     cout << "Digit Id " << " Amplitude " <<  " Index "  <<  " Nprim " << " Primaries list " <<  endl;      
1751
1752
1753     Int_t index ;
1754
1755
1756     for (index = 0 ; index < fDigits->GetEntries() ; index++) {
1757
1758
1759       digit = (AliEMCALDigit * )  fDigits->At(index) ;
1760
1761
1762       cout << setw(8)  <<  digit->GetId() << " "  <<    setw(3)  <<  digit->GetAmp() <<   "  "  
1763
1764
1765            << setw(6)  <<  digit->GetIndexInList() << "  "   
1766
1767
1768            << setw(5)  <<  digit->GetNprimary() <<"  ";
1769
1770
1771       
1772
1773
1774       Int_t iprimary;
1775
1776
1777       for (iprimary=0; iprimary<digit->GetNprimary(); iprimary++)
1778
1779
1780         cout << setw(5)  <<  digit->GetPrimary(iprimary+1) << " ";
1781
1782
1783       cout << endl;      
1784
1785
1786     }
1787
1788
1789     
1790
1791
1792   }
1793
1794
1795 }
1796
1797
1798 //_________________________________________________________________________________________
1799
1800
1801 void AliEMCALDigitizer::WriteDigits(Int_t event)
1802
1803
1804 {
1805
1806
1807
1808
1809
1810   // Makes TreeD in the output file. 
1811
1812
1813   // Check if branch already exists: 
1814
1815
1816   //   if yes, exit without writing: ROOT TTree does not support overwriting/updating of 
1817
1818
1819   //      already existing branches. 
1820
1821
1822   //   else creates branch with Digits, named "EMCAL", title "...",
1823
1824
1825   //      and branch "AliEMCALDigitizer", with the same title to keep all the parameters
1826
1827
1828   //      and names of files, from which digits are made.
1829
1830
1831
1832
1833
1834   AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; 
1835
1836
1837   const TClonesArray * digits = gime->Digits(GetName()) ; 
1838
1839
1840  TTree * treeD ;
1841
1842
1843
1844
1845
1846   if(fManager)
1847
1848
1849     treeD = fManager->GetTreeD() ;
1850
1851
1852  else
1853
1854
1855     treeD = gAlice->TreeD();
1856
1857
1858   
1859
1860
1861   // create new branches
1862
1863
1864   // -- generate file name if necessary
1865
1866
1867   char * file =0;
1868
1869
1870   if(gSystem->Getenv("CONFIG_SPLIT_FILE")){ //generating file name
1871
1872
1873     file = new char[strlen(gAlice->GetBaseFile())+20] ;
1874
1875
1876     sprintf(file,"%s/EMCAL.Digits.root",gAlice->GetBaseFile()) ;
1877
1878
1879   }
1880
1881
1882
1883
1884
1885   TDirectory *cwd = gDirectory;
1886
1887
1888   // -- create Digits branch
1889
1890
1891   Int_t bufferSize = 32000 ;    
1892
1893
1894   TBranch * digitsBranch = treeD->Branch("EMCAL",&digits,bufferSize);
1895
1896
1897   digitsBranch->SetTitle(GetName());
1898
1899
1900   if (file) {
1901
1902
1903     digitsBranch->SetFile(file);
1904
1905
1906     TIter next( digitsBranch->GetListOfBranches());
1907
1908
1909     TBranch * sbr ;
1910
1911
1912     while ((sbr=(TBranch*)next())) {
1913
1914
1915       sbr->SetFile(file);
1916
1917
1918     }   
1919
1920
1921     cwd->cd();
1922
1923
1924   } 
1925
1926
1927     
1928
1929
1930   // -- Create Digitizer branch
1931
1932
1933   Int_t splitlevel = 0 ;
1934
1935
1936   AliEMCALDigitizer * d = gime->Digitizer(GetName()) ;
1937
1938
1939   TBranch * digitizerBranch = treeD->Branch("AliEMCALDigitizer", "AliEMCALDigitizer", &d,bufferSize,splitlevel); 
1940
1941
1942   digitizerBranch->SetTitle(GetName());
1943
1944
1945   if (file) {
1946
1947
1948     digitizerBranch->SetFile(file);
1949
1950
1951     TIter next( digitizerBranch->GetListOfBranches());
1952
1953
1954     TBranch * sbr;
1955
1956
1957     while ((sbr=(TBranch*)next())) {
1958
1959
1960       sbr->SetFile(file);
1961
1962
1963     }   
1964
1965
1966     cwd->cd();
1967
1968
1969   }
1970
1971
1972   
1973
1974
1975   digitsBranch->Fill() ;      
1976
1977
1978   digitizerBranch->Fill() ;
1979
1980
1981
1982
1983
1984   treeD->Write(0,kOverwrite) ;  
1985
1986
1987  
1988
1989
1990 }
1991
1992
1993 //____________________________________________________________________________ 
1994
1995
1996 Float_t AliEMCALDigitizer::FrontEdgeTime(TClonesArray * ticks) 
1997
1998
1999 { // 
2000
2001
2002   ticks->Sort() ; //Sort in accordance with times of ticks
2003
2004
2005   TIter it(ticks) ;
2006
2007
2008   AliEMCALTick * ctick = (AliEMCALTick *) it.Next() ;
2009
2010
2011   Float_t time = ctick->CrossingTime(fTimeThreshold) ;    
2012
2013
2014   
2015
2016
2017   AliEMCALTick * t ;  
2018
2019
2020   while((t=(AliEMCALTick*) it.Next())){
2021
2022
2023     if(t->GetTime() < time)  //This tick starts before crossing
2024
2025
2026       *ctick+=*t ;
2027
2028
2029     else
2030
2031
2032       return time ;
2033
2034
2035     
2036
2037
2038     time = ctick->CrossingTime(fTimeThreshold) ;    
2039
2040
2041   }
2042
2043
2044   return time ;
2045
2046
2047 }
2048
2049
2050 //____________________________________________________________________________ 
2051
2052
2053 Float_t AliEMCALDigitizer::TimeOfNoise(void)
2054
2055
2056 {  // Calculates the time signal generated by noise
2057
2058
2059   //to be rewritten, now returns just big number
2060
2061
2062   return 1. ;
2063
2064
2065
2066
2067
2068 }
2069
2070
2071 //____________________________________________________________________________
2072
2073
2074 void AliEMCALDigitizer::SetSDigitsBranch(const char* title)
2075
2076
2077 {
2078
2079
2080   // we set title (comment) of the SDigits branch in the first! header file
2081
2082
2083   if( strcmp(GetName(), "") == 0 )
2084
2085
2086     Init() ;
2087
2088
2089
2090
2091
2092   AliEMCALGetter::GetInstance()->SDigits()->SetName(title) ; 
2093
2094
2095 }
2096
2097