]> git.uio.no Git - u/mrichter/AliRoot.git/blame - EMCAL/AliEMCALTrigger.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Warnings
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALTrigger.cxx
CommitLineData
f0377b23 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/* $Id$ */
f0377b23 16
17//_________________________________________________________________________
18//
19// Class for trigger analysis.
0964c2e9 20// Digits are grouped in TRU's (Trigger Units). A TRU consists of 384
85c25c2e 21// modules ordered fNTRUPhi x fNTRUEta. The algorithm searches all possible 2x2
0964c2e9 22// and nxn (n is a multiple of 2) cell combinations per each TRU, adding the
23// digits amplitude and finding the maximum. If found, look if it is isolated.
24// Maxima are transformed in ADC time samples. Each time bin is compared to the trigger
25// threshold until it is larger and then, triggers are set. Thresholds need to be fixed.
0b2ec9f7 26// Thresholds need to be fixed. Last 2 modules are half size in Phi, I considered
27// that the number of TRU is maintained for the last modules but decision not taken.
28// If different, then this must be changed.
f0377b23 29// Usage:
30//
31// //Inside the event loop
32// AliEMCALTrigger *tr = new AliEMCALTrigger();//Init Trigger
0b2ec9f7 33// tr->SetL0Threshold(100); //Arbitrary threshold values
85c25c2e 34// tr->SetL1GammaLowPtThreshold(1000);
35// tr->SetL1GammaMediumPtThreshold(10000);
36// tr->SetL1GammaHighPtThreshold(20000);
0964c2e9 37// ...
f0377b23 38// tr->Trigger(); //Execute Trigger
39// tr->Print(""); //Print results
40//
41//*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
42//////////////////////////////////////////////////////////////////////////////
43
d1f9a8ab 44#include <cassert>
f0377b23 45
46// --- ROOT system ---
85c25c2e 47#include <TTree.h>
48#include <TBranch.h>
49#include <TBrowser.h>
50#include <TH2F.h>
f0377b23 51
52// --- ALIROOT system ---
f0377b23 53#include "AliRun.h"
54#include "AliRunLoader.h"
55#include "AliTriggerInput.h"
56#include "AliEMCAL.h"
57#include "AliEMCALLoader.h"
58#include "AliEMCALDigit.h"
59#include "AliEMCALTrigger.h"
60#include "AliEMCALGeometry.h"
133abe1f 61#include "AliEMCALRawUtils.h"
0293e426 62#include "AliLog.h"
92d9f317 63#include "AliCaloConstants.h"
64
65using namespace CALO;
f0377b23 66
67ClassImp(AliEMCALTrigger)
68
85c25c2e 69TString AliEMCALTrigger::fgNameOfJetTriggers("EMCALJetTriggerL1");
70
f0377b23 71//______________________________________________________________________
72AliEMCALTrigger::AliEMCALTrigger()
9946f2fe 73 : AliTriggerDetector(), fGeom(0),
85c25c2e 74 f2x2MaxAmp(-1), f2x2ModulePhi(-1), f2x2ModuleEta(-1),
18a21c7c 75 f2x2SM(0),
85c25c2e 76 fnxnMaxAmp(-1), fnxnModulePhi(-1), fnxnModuleEta(-1),
0b2ec9f7 77 fnxnSM(0),
0964c2e9 78 fADCValuesHighnxn(0),fADCValuesLownxn(0),
79 fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
80 fDigitsList(0),
85c25c2e 81 fL0Threshold(100),fL1GammaLowPtThreshold(200),
82 fL1GammaMediumPtThreshold(500), fL1GammaHighPtThreshold(1000),
0964c2e9 83 fPatchSize(1), fIsolPatchSize(1),
84 f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1),
85 f2x2AmpOutOfPatchThres(100000), fnxnAmpOutOfPatchThres(100000),
86 fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),
85c25c2e 87 fSimulation(kTRUE), fIsolateInSuperModule(kTRUE), fTimeKey(kFALSE),
88 fAmpTrus(0),fTimeRtrus(0),fAmpSMods(0),
89 fTriggerPosition(6), fTriggerAmplitudes(4),
90 fNJetPatchPhi(3), fNJetPatchEta(3), fNJetThreshold(3), fL1JetThreshold(0), fJetMaxAmp(0),
91 fAmpJetMatrix(0), fJetMatrixE(0), fAmpJetMax(6,1), fVZER0Mult(0.)
f0377b23 92{
59264fa6 93 //ctor
0964c2e9 94 fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins];
95 fADCValuesLownxn = 0x0; //new Int_t[fTimeBins];
96 fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins];
97 fADCValuesLow2x2 = 0x0; //new Int_t[fTimeBins];
59264fa6 98
59264fa6 99 SetName("EMCAL");
85c25c2e 100 // Define jet threshold - can not change from outside now
101 // fNJetThreshold = 7; // For MB Pythia suppression
102 fNJetThreshold = 10; // Hijing
103 fL1JetThreshold = new Double_t[fNJetThreshold];
104 if(fNJetThreshold == 7) {
105 fL1JetThreshold[0] = 5./0.0153;
106 fL1JetThreshold[1] = 8./0.0153;
107 fL1JetThreshold[2] = 10./0.0153;
108 fL1JetThreshold[3] = 12./0.0153;
109 fL1JetThreshold[4] = 13./0.0153;
110 fL1JetThreshold[5] = 14./0.0153;
111 fL1JetThreshold[6] = 15./0.0153;
112 } else if(fNJetThreshold == 10) {
113 Double_t thGev[10]={5.,8.,10., 12., 13.,14.,15., 17., 20., 25.};
114 for(Int_t i=0; i<fNJetThreshold; i++) fL1JetThreshold[i] = thGev[i]/0.0153;
115 } else {
116 fL1JetThreshold[0] = 5./0.0153;
117 fL1JetThreshold[1] = 10./0.0153;
118 fL1JetThreshold[2] = 15./0.0153;
119 fL1JetThreshold[3] = 20./0.0153;
120 fL1JetThreshold[4] = 25./0.0153;
121 }
122 //
59264fa6 123 CreateInputs();
85c25c2e 124
125 fInputs.SetName("TriggersInputs");
59264fa6 126 //Print("") ;
f0377b23 127}
128
129
130
131//____________________________________________________________________________
132AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig)
18a21c7c 133 : AliTriggerDetector(trig),
9946f2fe 134 fGeom(trig.fGeom),
18a21c7c 135 f2x2MaxAmp(trig.f2x2MaxAmp),
85c25c2e 136 f2x2ModulePhi(trig.f2x2ModulePhi),
137 f2x2ModuleEta(trig.f2x2ModuleEta),
18a21c7c 138 f2x2SM(trig.f2x2SM),
0b2ec9f7 139 fnxnMaxAmp(trig.fnxnMaxAmp),
85c25c2e 140 fnxnModulePhi(trig.fnxnModulePhi),
141 fnxnModuleEta(trig.fnxnModuleEta),
0b2ec9f7 142 fnxnSM(trig.fnxnSM),
143 fADCValuesHighnxn(trig.fADCValuesHighnxn),
144 fADCValuesLownxn(trig.fADCValuesLownxn),
18a21c7c 145 fADCValuesHigh2x2(trig.fADCValuesHigh2x2),
146 fADCValuesLow2x2(trig.fADCValuesLow2x2),
147 fDigitsList(trig.fDigitsList),
148 fL0Threshold(trig.fL0Threshold),
85c25c2e 149 fL1GammaLowPtThreshold(trig.fL1GammaLowPtThreshold),
150 fL1GammaMediumPtThreshold(trig.fL1GammaMediumPtThreshold),
151 fL1GammaHighPtThreshold(trig.fL1GammaHighPtThreshold),
0964c2e9 152 fPatchSize(trig.fPatchSize),
153 fIsolPatchSize(trig.fIsolPatchSize),
154 f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch),
155 fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch),
156 f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres),
157 fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres),
158 fIs2x2Isol(trig.fIs2x2Isol),
159 fIsnxnIsol(trig.fIsnxnIsol),
160 fSimulation(trig.fSimulation),
85c25c2e 161 fIsolateInSuperModule(trig.fIsolateInSuperModule),
162 fTimeKey(trig.fTimeKey),
163 fAmpTrus(trig.fAmpTrus),
164 fTimeRtrus(trig.fTimeRtrus),
165 fAmpSMods(trig.fAmpSMods),
166 fTriggerPosition(trig.fTriggerPosition),
167 fTriggerAmplitudes(trig.fTriggerAmplitudes),
168 fNJetPatchPhi(trig.fNJetPatchPhi),
169 fNJetPatchEta(trig.fNJetPatchEta),
170 fNJetThreshold(trig.fNJetThreshold),
171 fL1JetThreshold(trig.fL1JetThreshold),
172 fJetMaxAmp(trig.fJetMaxAmp),
173 fAmpJetMatrix(trig.fAmpJetMatrix),
174 fJetMatrixE(trig.fJetMatrixE),
175 fAmpJetMax(trig.fAmpJetMax),
176 fVZER0Mult(trig.fVZER0Mult)
f0377b23 177{
f0377b23 178 // cpy ctor
f0377b23 179}
180
0293e426 181//____________________________________________________________________________
c35bbfd4 182AliEMCALTrigger::~AliEMCALTrigger() {
581ee158 183
184 //dtor
185
85c25c2e 186 if(GetTimeKey()) {
187 delete [] fADCValuesHighnxn;
188 delete [] fADCValuesLownxn;
189 delete [] fADCValuesHigh2x2;
190 delete [] fADCValuesLow2x2;
191 }
192 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
193 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
194 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
195 if(fAmpJetMatrix) delete fAmpJetMatrix;
196 if(fJetMatrixE) delete fJetMatrixE;
197 if(fL1JetThreshold) delete [] fL1JetThreshold;
c35bbfd4 198}
199
f0377b23 200//----------------------------------------------------------------------
201void AliEMCALTrigger::CreateInputs()
202{
203 // inputs
204
205 // Do not create inputs again!!
206 if( fInputs.GetEntriesFast() > 0 ) return;
85c25c2e 207
208 // Second parameter should be detector name = "EMCAL"
209 TString det("EMCAL"); // Apr 29, 2008
210 fInputs.AddLast( new AliTriggerInput( det+"_L0", det, 0x02) );
211 fInputs.AddLast( new AliTriggerInput( det+"_GammaHPt_L1", det, 0x04 ) );
212 fInputs.AddLast( new AliTriggerInput( det+"_GammaMPt_L1", det, 0x08 ) );
213 fInputs.AddLast( new AliTriggerInput( det+"_GammaLPt_L1", det, 0x016 ) );
ad534ff6 214 fInputs.AddLast( new AliTriggerInput( det+"_JetHPt_L1", det, 0x032 ) );
2994ff3b 215 fInputs.AddLast( new AliTriggerInput( det+"_JetMPt_L1", det, 0x048 ) );
216 fInputs.AddLast( new AliTriggerInput( det+"_JetLPt_L1", det, 0x064 ) );
85c25c2e 217
218 if(fNJetThreshold<=0) return;
219 // Jet Trigger(s)
220 UInt_t level = 0x032;
221 for(Int_t i=0; i<fNJetThreshold; i++ ) {
7e42cd80 222 TString name(GetNameOfJetTrigger(i));
85c25c2e 223 TString title("EMCAL Jet triger L1 :"); // unused now
224 // 0.0153 - hard coded now
225 title += Form("Th %i(%5.1f GeV) :", (Int_t)fL1JetThreshold[i], fL1JetThreshold[i] * 0.0153);
226 title += Form("patch %ix%i~(%3.2f(#phi)x%3.2f(#eta)) ",
227 fNJetPatchPhi, fNJetPatchEta, 0.11*(fNJetPatchPhi), 0.11*(fNJetPatchEta));
228 fInputs.AddLast( new AliTriggerInput(name, det, UChar_t(level)) );
229 level *= 2;
230 }
f0377b23 231
232}
233
234//____________________________________________________________________________
0964c2e9 235Bool_t AliEMCALTrigger::IsPatchIsolated(Int_t iPatchType, const TClonesArray * ampmatrixes, const Int_t iSM, const Int_t mtru, const Float_t maxamp, const Int_t maxphi, const Int_t maxeta) {
236
85c25c2e 237 // Nov 8, 2007
238 // EMCAL RTU size is 4modules(phi) x 24modules (eta)
239 // So maximum size of patch is 4modules x 4modules (EMCAL L0 trigger).
240 // Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch,
241 // inside isolation patch . iPatchType = 0 means calculation for 2x2 patch,
242 // iPatchType = 1 means calculation for nxn patch.
243 // In the next table there is an example of the different options of patch size and isolation patch size:
244 // Patch Size (fPatchSize)
245 // 0 1
246 // fIsolPatchSize 0 2x2 (not overlap) 4x4 (overlapped)
247 // 1 4x4 8x8
0964c2e9 248
249 Bool_t b = kFALSE;
622e10be 250 if(!ampmatrixes) return kFALSE;
251
85c25c2e 252 // Get matrix of TRU or Module with maximum amplitude patch.
253 Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36.
0964c2e9 254 TMatrixD * ampmatrix = 0x0;
255 Int_t colborder = 0;
256 Int_t rowborder = 0;
85c25c2e 257 static int keyPrint = 0;
0293e426 258 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated : iSM %i mtru %i itru %i maxphi %i maxeta %i \n", iSM, mtru, itru, maxphi, maxeta));
0964c2e9 259
85c25c2e 260 if(fIsolateInSuperModule){ // ?
0964c2e9 261 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
85c25c2e 262 rowborder = fGeom->GetNPhi();
263 colborder = fGeom->GetNZ();
0964c2e9 264 AliDebug(2,"Isolate trigger in Module");
85c25c2e 265 } else{
0964c2e9 266 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
85c25c2e 267 rowborder = fGeom->GetNModulesInTRUPhi();
268 colborder = fGeom->GetNModulesInTRUEta();
0964c2e9 269 AliDebug(2,"Isolate trigger in TRU");
270 }
85c25c2e 271 if(iSM>9) rowborder /= 2; // half size in phi
0964c2e9 272
6ce4e2ed 273 if(!ampmatrixes || !ampmatrix){
a51e676d 274 AliError("Could not recover the matrix with the amplitudes");
275 return kFALSE;
276 }
277
85c25c2e 278 //Define patch modules - what is this ??
279 Int_t isolmodules = fIsolPatchSize*(1+iPatchType);
280 Int_t ipatchmodules = 2*(1+fPatchSize*iPatchType);
281 Int_t minrow = maxphi - isolmodules;
282 Int_t mincol = maxeta - isolmodules;
283 Int_t maxrow = maxphi + isolmodules + ipatchmodules;
284 Int_t maxcol = maxeta + isolmodules + ipatchmodules;
285
286 minrow = minrow<0?0 :minrow;
287 mincol = mincol<0?0 :mincol;
288
289 maxrow = maxrow>rowborder?rowborder :maxrow;
290 maxcol = maxcol>colborder?colborder :maxcol;
9946f2fe 291
85c25c2e 292 //printf("%s\n",Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
293 //printf("%s\n",Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
294 // AliDebug(2,Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
295 //AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
0964c2e9 296
0964c2e9 297 //Add amplitudes in all isolation patch
85c25c2e 298 Float_t amp = 0.;
0964c2e9 299 for(Int_t irow = minrow ; irow < maxrow; irow ++)
300 for(Int_t icol = mincol ; icol < maxcol ; icol ++)
301 amp += (*ampmatrix)(irow,icol);
302
303 AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
304
305 if(amp < maxamp){
85c25c2e 306 // AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
307 // ampmatrix->Print();
0964c2e9 308 return kFALSE;
85c25c2e 309 } else {
0964c2e9 310 amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
85c25c2e 311 }
0964c2e9 312
313 AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
314
315 //Fill isolation amplitude data member and say if patch is isolated.
316 if(iPatchType == 0){ //2x2 case
317 f2x2AmpOutOfPatch = amp;
85c25c2e 318 if(amp < f2x2AmpOutOfPatchThres) b=kTRUE;
319 } else if(iPatchType == 1){ //nxn case
0964c2e9 320 fnxnAmpOutOfPatch = amp;
85c25c2e 321 if(amp < fnxnAmpOutOfPatchThres) b=kTRUE;
0964c2e9 322 }
323
0293e426 324 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated - OUT \n"));
85c25c2e 325
0964c2e9 326 return b;
327
328}
329
85c25c2e 330/*
0964c2e9 331//____________________________________________________________________________
c35bbfd4 332void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD &ampmax2, TMatrixD &ampmaxn){
f0377b23 333
85c25c2e 334 //Sums energy of all possible 2x2 (L0) and nxn (L1) modules per each TRU.
335 //Fast signal in the experiment is given by 2x2 modules,
336 //for this reason we loop inside the TRU modules by 2.
33d0b833 337
59264fa6 338 //Declare and initialize variables
9946f2fe 339 Int_t nCellsPhi = fGeom->GetNCellsInTRUPhi();
33d0b833 340 if(isupermod > 9)
59264fa6 341 nCellsPhi = nCellsPhi / 2 ; //Half size SM. Not Final.
f0377b23 342 // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
9946f2fe 343 Int_t nCellsEta = fGeom->GetNCellsInTRUEta();
344 Int_t nTRU = fGeom->GetNTRU();
f0377b23 345 // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
0964c2e9 346 //Int_t nTRU = geom->GeNTRU();//3 TRU per super module
f0377b23 347
59264fa6 348 Float_t amp2 = 0 ;
0b2ec9f7 349 Float_t ampn = 0 ;
33d0b833 350 for(Int_t i = 0; i < 4; i++){
9946f2fe 351 for(Int_t j = 0; j < nTRU; j++){
c35bbfd4 352 ampmax2(i,j) = -1;
353 ampmaxn(i,j) = -1;
59264fa6 354 }
355 }
f0377b23 356
59264fa6 357 //Create matrix that will contain 2x2 amplitude sums
0b2ec9f7 358 //used to calculate the nxn sums
c35bbfd4 359 TMatrixD tru2x2(nCellsPhi/2,nCellsEta/2) ;
59264fa6 360 for(Int_t i = 0; i < nCellsPhi/2; i++)
361 for(Int_t j = 0; j < nCellsEta/2; j++)
c35bbfd4 362 tru2x2(i,j) = -1;
59264fa6 363
364 //Loop over all TRUS in a supermodule
9946f2fe 365 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
59264fa6 366 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
367 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
9946f2fe 368 Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
33d0b833 369
59264fa6 370 //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
371 for(Int_t irow = 0 ; irow < nCellsPhi; irow += 2){
372 for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
373 amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
374 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
33d0b833 375
0964c2e9 376 //Fill matrix with added 2x2 cells for use in nxn sums
c35bbfd4 377 tru2x2(irow/2,icol/2) = amp2 ;
59264fa6 378 //Select 2x2 maximum sums to select L0
c35bbfd4 379 if(amp2 > ampmax2(0,mtru)){
380 ampmax2(0,mtru) = amp2 ;
381 ampmax2(1,mtru) = irow;
382 ampmax2(2,mtru) = icol;
59264fa6 383 }
384 }
385 }
386
387 //Find most recent time in the selected 2x2 cell
c35bbfd4 388 ampmax2(3,mtru) = 1 ;
389 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
390 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
59264fa6 391 for(Int_t i = 0; i<2; i++){
392 for(Int_t j = 0; j<2; j++){
393 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
c35bbfd4 394 if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) )
395 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j);
59264fa6 396 }
397 }
398 }
0b2ec9f7 399
400 //Sliding nxn, add nxn amplitudes (OVERLAP)
401 if(fPatchSize > 0){
402 for(Int_t irow = 0 ; irow < nCellsPhi/2; irow++){
403 for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
404 ampn = 0;
405 if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
406 for(Int_t i = 0 ; i <= fPatchSize ; i++)
407 for(Int_t j = 0 ; j <= fPatchSize ; j++)
c35bbfd4 408 ampn += tru2x2(irow+i,icol+j);
0b2ec9f7 409 //Select nxn maximum sums to select L1
c35bbfd4 410 if(ampn > ampmaxn(0,mtru)){
411 ampmaxn(0,mtru) = ampn ;
412 ampmaxn(1,mtru) = irow*2;
413 ampmaxn(2,mtru) = icol*2;
0b2ec9f7 414 }
59264fa6 415 }
416 }
417 }
0b2ec9f7 418
419 //Find most recent time in selected nxn cell
c35bbfd4 420 ampmaxn(3,mtru) = 1 ;
421 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
422 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
0b2ec9f7 423 for(Int_t i = 0; i<4*fPatchSize; i++){
424 for(Int_t j = 0; j<4*fPatchSize; j++){
9946f2fe 425 if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta){//Avoid exit the TRU
0b2ec9f7 426 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
c35bbfd4 427 if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) )
428 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j);
0b2ec9f7 429 }
430 }
59264fa6 431 }
432 }
433 }
0b2ec9f7 434 else {
c35bbfd4 435 ampmaxn(0,mtru) = ampmax2(0,mtru);
436 ampmaxn(1,mtru) = ampmax2(1,mtru);
437 ampmaxn(2,mtru) = ampmax2(2,mtru);
438 ampmaxn(3,mtru) = ampmax2(3,mtru);
0b2ec9f7 439 }
f0377b23 440 }
f0377b23 441}
85c25c2e 442*/
443//____________________________________________________________________________
444void AliEMCALTrigger::MakeSlidingTowers(const TClonesArray * amptrus, const TClonesArray * timeRtrus,
a51e676d 445 const Int_t isupermod,TMatrixD &ampmax2, TMatrixD &ampmaxn){
85c25c2e 446
447 // Output from module (2x2 cells from one module)
448 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi(); // now 4 modules (3 div in phi)
449 if(isupermod > 9)
450 nModulesPhi = nModulesPhi / 2 ; // Half size SM. Not Final.
451 //
452 Int_t nModulesEta = fGeom->GetNModulesInTRUEta(); // now 24 modules (no division in eta)
453 Int_t nTRU = fGeom->GetNTRU();
454 static int keyPrint = 0;
0293e426 455 if(keyPrint) AliDebug(2,Form("MakeSlidingTowers : nTRU %i nModulesPhi %i nModulesEta %i ",
a51e676d 456 nTRU, nModulesPhi, nModulesEta ));
457
85c25c2e 458 Float_t amp2 = 0 ;
459 Float_t ampn = 0 ;
460 for(Int_t i = 0; i < 4; i++){
461 for(Int_t j = 0; j < nTRU; j++){
462 ampmax2(i,j) = ampmaxn(i,j) = -1;
463 }
464 }
a51e676d 465
85c25c2e 466 // Create matrix that will contain 2x2 amplitude sums
467 // used to calculate the nxn sums
468 TMatrixD tru2x2(nModulesPhi/2,nModulesEta/2);
a51e676d 469
85c25c2e 470 // Loop over all TRUS in a supermodule
471 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
472 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
473 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
474 Int_t mtru = itru - isupermod*nTRU ; // Number of TRU in Supermodule !!
a51e676d 475
476 if(!amptru || !timeRtru){
477 AliError("Amplitude or Time TRU matrix not available")
478 return;
479 }
480
85c25c2e 481 // Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
482 for(Int_t irow = 0 ; irow < nModulesPhi; irow +=2){
483 for(Int_t icol = 0 ; icol < nModulesEta ; icol +=2){
a51e676d 484 amp2 = (*amptru)(irow,icol) +(*amptru)(irow+1,icol)+
485 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
486
487 //Fill matrix with added 2x2 towers for use in nxn sums
488 tru2x2(irow/2,icol/2) = amp2 ;
489 //Select 2x2 maximum sums to select L0
490 if(amp2 > ampmax2(0,mtru)){
491 ampmax2(0,mtru) = amp2 ;
492 ampmax2(1,mtru) = irow;
493 ampmax2(2,mtru) = icol;
494 }
85c25c2e 495 }
496 }
497
498 ampmax2(3,mtru) = 0.;
499 if(GetTimeKey()) {
a51e676d 500 // Find most recent time in the selected 2x2 towers
85c25c2e 501 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
502 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
503 for(Int_t i = 0; i<2; i++){
504 for(Int_t j = 0; j<2; j++){
a51e676d 505 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
506 if((*timeRtru)(row2+i,col2+j) > ampmax2(3,mtru) )
507 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); // max time
508 }
85c25c2e 509 }
510 }
511 }
a51e676d 512
85c25c2e 513 //Sliding nxn, add nxn amplitudes (OVERLAP)
514 if(fPatchSize > 0){
515 for(Int_t irow = 0 ; irow < nModulesPhi/2; irow++){
a51e676d 516 for(Int_t icol = 0 ; icol < nModulesEta/2; icol++){
517 ampn = 0;
518 if( (irow+fPatchSize) < nModulesPhi/2 && (icol+fPatchSize) < nModulesEta/2){ //Avoid exit the TRU
519 for(Int_t i = 0 ; i <= fPatchSize ; i++)
520 for(Int_t j = 0 ; j <= fPatchSize ; j++)
521 ampn += tru2x2(irow+i,icol+j);
522 //Select nxn maximum sums to select L1
523 if(ampn > ampmaxn(0,mtru)){
524 ampmaxn(0,mtru) = ampn ;
525 ampmaxn(1,mtru) = irow;
526 ampmaxn(2,mtru) = icol;
527 }
528 }
529 }
85c25c2e 530 }
531
532 ampmaxn(3,mtru) = 0.; // Was 1 , I don't know why
533 if(GetTimeKey()) {
a51e676d 534 //Find most recent time in selected nxn cell
85c25c2e 535 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
536 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
537 for(Int_t i = 0; i<4*fPatchSize; i++){
a51e676d 538 for(Int_t j = 0; j<4*fPatchSize; j++){
539 if( (rown+i) < nModulesPhi && (coln+j) < nModulesEta){//Avoid exit the TRU
540 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
541 if((*timeRtru)(rown+i,coln+j) > ampmaxn(3,mtru) )
542 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); // max time
543 }
544 }
545 }
85c25c2e 546 }
547 }
548 } else { // copy 2x2 to nxn
a51e676d 549 ampmaxn(0,mtru) = ampmax2(0,mtru);
550 ampmaxn(1,mtru) = ampmax2(1,mtru);
551 ampmaxn(2,mtru) = ampmax2(2,mtru);
552 ampmaxn(3,mtru) = ampmax2(3,mtru);
85c25c2e 553 }
554 }
0293e426 555 if(keyPrint) AliDebug(2,Form(" : MakeSlidingTowers -OUt \n"));
85c25c2e 556}
f0377b23 557
558//____________________________________________________________________________
559void AliEMCALTrigger::Print(const Option_t * opt) const
560{
561
562 //Prints main parameters
563
564 if(! opt)
565 return;
566 AliTriggerInput* in = 0x0 ;
0293e426 567 AliInfo(Form(" fSimulation %i (input option) : #digits %i\n", fSimulation, fDigitsList->GetEntries()));
568 AliInfo(Form(" fTimeKey %i \n ", fTimeKey));
569
570 AliInfo(Form("\t Maximum Amplitude after Sliding Cell, \n")) ;
571 AliInfo(Form("\t -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
572 f2x2MaxAmp,f2x2SM)) ;
573 AliInfo(Form("\t -2x2 from row %d to row %d and from column %d to column %d\n", f2x2ModulePhi, f2x2ModulePhi+2, f2x2ModuleEta, f2x2ModuleEta+2));
574 AliInfo(Form("\t -2x2 Isolation Patch %d x %d, Amplitude out of 2x2 patch is %f, threshold %f, Isolated? %d \n", 2*fIsolPatchSize+2, 2*fIsolPatchSize+2, f2x2AmpOutOfPatch, f2x2AmpOutOfPatchThres,static_cast<Int_t> (fIs2x2Isol)));
0b2ec9f7 575 if(fPatchSize > 0){
0293e426 576 AliInfo(Form("\t Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1)));
577 AliInfo(Form("\t -nxn cells sum (overlapped) : %10.2f, in Super Module %d\n", fnxnMaxAmp,fnxnSM));
578 AliInfo(Form("\t -nxn from row %d to row %d and from column %d to column %d\n", fnxnModulePhi, fnxnModulePhi+4*fPatchSize, fnxnModuleEta, fnxnModuleEta+4*fPatchSize)) ;
579 AliInfo(Form("\t -nxn Isolation Patch %d x %d, Amplitude out of nxn patch is %f, threshold %f, Isolated? %d \n", 4*fIsolPatchSize+2*(fPatchSize+1),4*fIsolPatchSize+2*(fPatchSize+1) , fnxnAmpOutOfPatch, fnxnAmpOutOfPatchThres,static_cast<Int_t> (fIsnxnIsol) ));
0b2ec9f7 580 }
0964c2e9 581
0293e426 582 AliInfo(Form("\t Isolate in SuperModule? %d\n", fIsolateInSuperModule)) ;
583 AliInfo(Form("\t Threshold for LO %10.2f\n", fL0Threshold));
0964c2e9 584
59264fa6 585 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
f0377b23 586 if(in->GetValue())
0293e426 587 AliInfo(Form("\t *** EMCAL LO is set ***\n"));
f0377b23 588
0293e426 589 AliInfo(Form("\t Gamma Low Pt Threshold for L1 %10.2f\n", fL1GammaLowPtThreshold));
85c25c2e 590 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_GammaLPt_L1" );
f0377b23 591 if(in->GetValue())
0293e426 592 AliInfo(Form("\t *** EMCAL Gamma Low Pt for L1 is set ***\n"));
f0377b23 593
0293e426 594 AliInfo(Form("\t Gamma Medium Pt Threshold for L1 %10.2f\n", fL1GammaMediumPtThreshold));
85c25c2e 595 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaMPt_L1" );
f0377b23 596 if(in->GetValue())
0293e426 597 AliInfo(Form("\t *** EMCAL Gamma Medium Pt for L1 is set ***\n"));
f0377b23 598
0293e426 599 AliInfo(Form("\t Gamma High Pt Threshold for L1 %10.2f\n", fL1GammaHighPtThreshold));
85c25c2e 600 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaHPt_L1" );
f0377b23 601 if(in->GetValue())
0293e426 602 AliInfo(Form("\t *** EMCAL Gamma High Pt for L1 is set ***\n")) ;
f0377b23 603
604}
605
606//____________________________________________________________________________
0964c2e9 607void AliEMCALTrigger::SetTriggers(const TClonesArray * ampmatrix,const Int_t iSM,
c35bbfd4 608 const TMatrixD &ampmax2,
9946f2fe 609 const TMatrixD &ampmaxn)
f0377b23 610{
0b2ec9f7 611 //Checks the 2x2 and nxn maximum amplitude per each TRU and
59264fa6 612 //compares with the different L0 and L1 triggers thresholds
613 Float_t max2[] = {-1,-1,-1,-1} ;
0b2ec9f7 614 Float_t maxn[] = {-1,-1,-1,-1} ;
0964c2e9 615 Int_t mtru2 = -1 ;
616 Int_t mtrun = -1 ;
f0377b23 617
9946f2fe 618 Int_t nTRU = fGeom->GetNTRU();
619
59264fa6 620 //Find maximum summed amplitude of all the TRU
621 //in a Super Module
9946f2fe 622 for(Int_t i = 0 ; i < nTRU ; i++){
c35bbfd4 623 if(max2[0] < ampmax2(0,i) ){
624 max2[0] = ampmax2(0,i) ; // 2x2 summed max amplitude
625 max2[1] = ampmax2(1,i) ; // corresponding phi position in TRU
626 max2[2] = ampmax2(2,i) ; // corresponding eta position in TRU
627 max2[3] = ampmax2(3,i) ; // corresponding most recent time
0964c2e9 628 mtru2 = i ;
59264fa6 629 }
c35bbfd4 630 if(maxn[0] < ampmaxn(0,i) ){
631 maxn[0] = ampmaxn(0,i) ; // nxn summed max amplitude
632 maxn[1] = ampmaxn(1,i) ; // corresponding phi position in TRU
633 maxn[2] = ampmaxn(2,i) ; // corresponding eta position in TRU
634 maxn[3] = ampmaxn(3,i) ; // corresponding most recent time
0964c2e9 635 mtrun = i ;
59264fa6 636 }
637 }
638
639 //--------Set max amplitude if larger than in other Super Modules------------
640 Float_t maxtimeR2 = -1 ;
0b2ec9f7 641 Float_t maxtimeRn = -1 ;
133abe1f 642 static AliEMCALRawUtils rawUtil;
92d9f317 643 // Int_t nTimeBins = rawUtil.GetRawFormatTimeBins() ;
644 Int_t nTimeBins = TIMEBINS; //changed by PTH
59264fa6 645
646 //Set max of 2x2 amplitudes and select L0 trigger
647 if(max2[0] > f2x2MaxAmp ){
85c25c2e 648 // if(max2[0] > 5) printf(" L0 : iSM %i: max2[0] %5.0f : max2[3] %5.0f (maxtimeR2) \n",
649 // iSM, max2[0], max2[3]);
59264fa6 650 f2x2MaxAmp = max2[0] ;
651 f2x2SM = iSM ;
652 maxtimeR2 = max2[3] ;
85c25c2e 653 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtru2,
59264fa6 654 static_cast<Int_t>(max2[1]),
655 static_cast<Int_t>(max2[2]),
85c25c2e 656 f2x2ModulePhi,f2x2ModuleEta);
59264fa6 657
0964c2e9 658 //Isolated patch?
659 if(fIsolateInSuperModule)
85c25c2e 660 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, f2x2ModulePhi,f2x2ModuleEta) ;
0964c2e9 661 else
662 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, static_cast<Int_t>(max2[1]), static_cast<Int_t>(max2[2])) ;
663
85c25c2e 664 if(GetTimeKey()) {
59264fa6 665 //Transform digit amplitude in Raw Samples
85c25c2e 666 if (fADCValuesLow2x2 == 0) {
667 fADCValuesLow2x2 = new Int_t[nTimeBins];
668 fADCValuesHigh2x2 = new Int_t[nTimeBins];
669 }
670 //printf(" maxtimeR2 %12.5e (1)\n", maxtimeR2);
92d9f317 671 // rawUtil.RawSampledResponse(maxtimeR2 * AliEMCALRawUtils::GetRawFormatTimeBin(),
672 // f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
673
674 rawUtil.RawSampledResponse(maxtimeR2*TIMEBINMAX/TIMEBINS,
675 f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
676
85c25c2e 677 // Set Trigger Inputs, compare ADC time bins until threshold is attained
678 // Set L0
679 for(Int_t i = 0 ; i < nTimeBins ; i++){
680 // printf(" fADCValuesHigh2x2[%i] %i : %i \n", i, fADCValuesHigh2x2[i], fADCValuesLow2x2[i]);
681 if(fADCValuesHigh2x2[i] >= fL0Threshold || fADCValuesLow2x2[i] >= fL0Threshold){
682 SetInput("EMCAL_L0") ;
683 break;
684 }
685 }
686 } else {
687 // Nov 5 - no analysis of time information
688 if(f2x2MaxAmp >= fL0Threshold) { // should add the low amp too
689 SetInput("EMCAL_L0");
59264fa6 690 }
691 }
f0377b23 692 }
59264fa6 693
0b2ec9f7 694 //------------Set max of nxn amplitudes and select L1 trigger---------
695 if(maxn[0] > fnxnMaxAmp ){
696 fnxnMaxAmp = maxn[0] ;
697 fnxnSM = iSM ;
698 maxtimeRn = maxn[3] ;
85c25c2e 699 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtrun,
0b2ec9f7 700 static_cast<Int_t>(maxn[1]),
701 static_cast<Int_t>(maxn[2]),
85c25c2e 702 fnxnModulePhi,fnxnModuleEta) ;
0964c2e9 703
704 //Isolated patch?
705 if(fIsolateInSuperModule)
85c25c2e 706 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, fnxnModulePhi, fnxnModuleEta) ;
0964c2e9 707 else
708 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, static_cast<Int_t>(maxn[1]), static_cast<Int_t>(maxn[2])) ;
709
85c25c2e 710 if(GetTimeKey()) {
59264fa6 711 //Transform digit amplitude in Raw Samples
85c25c2e 712 if (fADCValuesLownxn == 0) {
713 fADCValuesHighnxn = new Int_t[nTimeBins];
714 fADCValuesLownxn = new Int_t[nTimeBins];
715 }
92d9f317 716 // rawUtil.RawSampledResponse(maxtimeRn * AliEMCALRawUtils::GetRawFormatTimeBin(),
717 // fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
718 rawUtil.RawSampledResponse(maxtimeRn*TIMEBINMAX/TIMEBINS,
85c25c2e 719 fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
92d9f317 720
721
59264fa6 722 //Set Trigger Inputs, compare ADC time bins until threshold is attained
723 //SetL1 Low
85c25c2e 724 for(Int_t i = 0 ; i < nTimeBins ; i++){
725 if(fADCValuesHighnxn[i] >= fL1GammaLowPtThreshold || fADCValuesLownxn[i] >= fL1GammaLowPtThreshold){
726 SetInput("EMCAL_GammaLPt_L1") ;
727 break;
728 }
59264fa6 729 }
59264fa6 730
731 //SetL1 Medium
85c25c2e 732 for(Int_t i = 0 ; i < nTimeBins ; i++){
733 if(fADCValuesHighnxn[i] >= fL1GammaMediumPtThreshold || fADCValuesLownxn[i] >= fL1GammaMediumPtThreshold){
734 SetInput("EMCAL_GammaMPt_L1") ;
735 break;
736 }
59264fa6 737 }
59264fa6 738
739 //SetL1 High
85c25c2e 740 for(Int_t i = 0 ; i < nTimeBins ; i++){
741 if(fADCValuesHighnxn[i] >= fL1GammaHighPtThreshold || fADCValuesLownxn[i] >= fL1GammaHighPtThreshold){
742 SetInput("EMCAL_GammaHPt_L1") ;
743 break;
744 }
745 }
746 } else {
747 // Nov 5 - no analysis of time information
748 if(fnxnMaxAmp >= fL1GammaLowPtThreshold) { // should add the low amp too
749 SetInput("EMCAL_GammaLPt_L1") ; //SetL1 Low
750 }
751 if(fnxnMaxAmp >= fL1GammaMediumPtThreshold) { // should add the low amp too
752 SetInput("EMCAL_GammaMPt_L1") ; //SetL1 Medium
753 }
754 if(fnxnMaxAmp >= fL1GammaHighPtThreshold) { // should add the low amp too
755 SetInput("EMCAL_GammaHPt_L1") ; //SetL1 High
59264fa6 756 }
757 }
85c25c2e 758 }
59264fa6 759}
f0377b23 760
c35bbfd4 761//____________________________________________________________________________
9946f2fe 762void AliEMCALTrigger::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * ampmatrixsmod, TClonesArray * timeRmatrix) {
c35bbfd4 763
764// Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
765// Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
766// TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
767// Last 2 modules are half size in Phi, I considered that the number of TRU
768// is maintained for the last modules but decision not taken. If different,
769// then this must be changed. Also fill a matrix with all amplitudes in supermodule for isolation studies.
770
85c25c2e 771// Initilize and declare variables
772// List of TRU matrices initialized to 0.
773// printf("<I> AliEMCALTrigger::FillTRU() started : # digits %i\n", digits->GetEntriesFast());
774
775// Nov 2, 2007.
776// One input per EMCAL module so size of matrix is reduced by 4 (2x2 division case)
777
778 Int_t nPhi = fGeom->GetNPhi();
779 Int_t nZ = fGeom->GetNZ();
780 Int_t nTRU = fGeom->GetNTRU();
781 // Int_t nTRUPhi = fGeom->GetNTRUPhi();
782 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi();
783 Int_t nModulesPhi2 = fGeom->GetNModulesInTRUPhi();
784 Int_t nModulesEta = fGeom->GetNModulesInTRUEta();
785 // printf("<I> AliEMCALTrigger::FillTRU() nTRU %i nTRUPhi %i : nModulesPhi %i nModulesEta %i \n",
786 // nTRU, nTRUPhi, nModulesPhi, nModulesEta);
c35bbfd4 787
788 Int_t id = -1;
789 Float_t amp = -1;
790 Float_t timeR = -1;
791 Int_t iSupMod = -1;
792 Int_t nModule = -1;
793 Int_t nIphi = -1;
794 Int_t nIeta = -1;
795 Int_t iphi = -1;
796 Int_t ieta = -1;
85c25c2e 797 // iphim, ietam - module indexes in SM
798 Int_t iphim = -1;
799 Int_t ietam = -1;
c35bbfd4 800
801 //List of TRU matrices initialized to 0.
9946f2fe 802 Int_t nSup = fGeom->GetNumberOfSuperModules();
803 for(Int_t k = 0; k < nTRU*nSup; k++){
85c25c2e 804 TMatrixD amptrus(nModulesPhi,nModulesEta) ;
805 TMatrixD timeRtrus(nModulesPhi,nModulesEta) ;
c35bbfd4 806 // Do we need to initialise? I think TMatrixD does it by itself...
85c25c2e 807 for(Int_t i = 0; i < nModulesPhi; i++){
808 for(Int_t j = 0; j < nModulesEta; j++){
c35bbfd4 809 amptrus(i,j) = 0.0;
810 timeRtrus(i,j) = 0.0;
811 }
812 }
813 new((*ampmatrix)[k]) TMatrixD(amptrus) ;
814 new((*timeRmatrix)[k]) TMatrixD(timeRtrus) ;
815 }
816
85c25c2e 817 // List of Modules matrices initialized to 0.
c35bbfd4 818 for(Int_t k = 0; k < nSup ; k++){
85c25c2e 819 int mphi = nPhi;
820 // if(nSup>9) mphi = nPhi/2; // the same size
821 TMatrixD ampsmods( mphi, nZ);
822 for(Int_t i = 0; i < mphi; i++){
823 for(Int_t j = 0; j < nZ; j++){
c35bbfd4 824 ampsmods(i,j) = 0.0;
825 }
826 }
85c25c2e 827 new((*ampmatrixsmod)[k]) TMatrixD(ampsmods) ;
c35bbfd4 828 }
829
830 AliEMCALDigit * dig ;
831
832 //Digits loop to fill TRU matrices with amplitudes.
833 for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
834
835 dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
a51e676d 836 if(dig){
837 amp = Float_t(dig->GetAmplitude()); // Energy of the digit (arbitrary units)
838 id = dig->GetId() ; // Id label of the cell
839 timeR = dig->GetTimeR() ; // Earliest time of the digit
840 if(amp<=0.0) AliDebug(1,Form(" id %i amp %f \n", id, amp));
841 // printf(" FILLTRU : timeR %10.5e time %10.5e : amp %10.5e \n", timeR, dig->GetTime(), amp);
842 // Get eta and phi cell position in supermodule
843 Bool_t bCell = fGeom->GetCellIndex(id, iSupMod, nModule, nIphi, nIeta) ;
844 if(!bCell)
845 AliError(Form("%i Wrong cell id number %i ", idig, id)) ;
846
847 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
848 // iphim, ietam - module indexes in SM
849 fGeom->GetModuleIndexesFromCellIndexesInSModule(iSupMod,iphi,ieta, iphim, ietam, nModule);
850 //if(iSupMod >9)
851 //printf("iSupMod %i nModule %i iphi %i ieta %i iphim %i ietam %i \n",
852 //iSupMod,nModule, iphi, ieta, iphim, ietam);
853
854 // Check to which TRU in the supermodule belongs the cell.
855 // Supermodules are divided in a TRU matrix of dimension
856 // (fNTRUPhi,fNTRUEta).
857 // Each TRU is a cell matrix of dimension (nModulesPhi,nModulesEta)
858
859 // First calculate the row and column in the supermodule
860 // of the TRU to which the cell belongs.
861 Int_t row = iphim / nModulesPhi;
862 Int_t col = ietam / nModulesEta;
863 //Calculate label number of the TRU
864 Int_t itru = fGeom->GetAbsTRUNumberFromNumberInSm(row, col, iSupMod);
865
866 //Fill TRU matrix with cell values
867 TMatrixD * amptrus = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
868 TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
869
6ce4e2ed 870 if(!amptrus || !timeRtrus){
a51e676d 871 AliError("Could not recover the TRU matrix with amplitudes or times");
a51e676d 872 }
622e10be 873 else{
874 //Calculate row and column of the module inside the TRU with number itru
875 Int_t irow = iphim - row * nModulesPhi;
876 if(iSupMod > 9)
877 irow = iphim - row * nModulesPhi2; // size of matrix the same
878 Int_t icol = ietam - col * nModulesEta;
a51e676d 879
622e10be 880 (*amptrus)(irow,icol) += amp ;
881 if((*timeRtrus)(irow,icol) <0.0 || (*timeRtrus)(irow,icol) <= timeR){ // ??
882 (*timeRtrus)(irow,icol) = timeR ;
883 }
a51e676d 884 }
885 //printf(" ieta %i iphi %i iSM %i || col %i row %i : itru %i -> amp %f\n",
886 // ieta, iphi, iSupMod, col, row, itru, amp);
887 //####################SUPERMODULE MATRIX ##################
888 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSupMod)) ;
889 if(!ampsmods){
890 AliError("Could not recover the matrix per SM");
891 continue;
892 }
893 (*ampsmods)(iphim,ietam) += amp ;
894 // printf(" id %i iphim %i ietam %i SM %i : irow %i icol %i itru %i : amp %6.0f\n",
895 //id, iphim, ietam, iSupMod, irow, icol, itru, amp);
85c25c2e 896 }
a51e676d 897 else AliError("Could not recover the digit");
c35bbfd4 898 }
85c25c2e 899 //assert(0);
900 //printf("<I> AliEMCALTrigger::FillTRU() is ended \n");
c35bbfd4 901}
0293e426 902
59264fa6 903//____________________________________________________________________________
904void AliEMCALTrigger::Trigger()
905{
906 //Main Method to select triggers.
581ee158 907 TH1::AddDirectory(0);
908
33c3c91a 909 AliRunLoader *runLoader = AliRunLoader::Instance();
85c25c2e 910 AliEMCALLoader *emcalLoader = 0;
911 if(runLoader) {
912 emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
913 }
0964c2e9 914
59264fa6 915 //Load EMCAL Geometry
a51e676d 916 if (runLoader && runLoader->GetAliRun()){
917 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"));
918 if(emcal)fGeom = emcal->GetGeometry();
919 }
920
921 if (!fGeom)
384c0bba 922 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
923
a51e676d 924 if (!fGeom)
59264fa6 925 AliFatal("Did not get geometry from EMCALLoader");
0964c2e9 926
59264fa6 927 //Define parameters
9946f2fe 928 Int_t nSuperModules = fGeom->GetNumberOfSuperModules() ; //12 SM in EMCAL
85c25c2e 929 Int_t nTRU = fGeom->GetNTRU(); // 3 TRU per super module
59264fa6 930
931 //Intialize data members each time the trigger is called in event loop
85c25c2e 932 f2x2MaxAmp = -1; f2x2ModulePhi = -1; f2x2ModuleEta = -1;
933 fnxnMaxAmp = -1; fnxnModulePhi = -1; fnxnModuleEta = -1;
59264fa6 934
85c25c2e 935 // Take the digits list if simulation
936 if(fSimulation && runLoader && emcalLoader){ // works than run seperate macros
c787fb51 937 runLoader->LoadDigits("EMCAL");
59264fa6 938 fDigitsList = emcalLoader->Digits() ;
85c25c2e 939 runLoader->LoadSDigits("EMCAL");
59264fa6 940 }
85c25c2e 941 // Digits list should be set by method SetDigitsList(TClonesArray * digits)
59264fa6 942 if(!fDigitsList)
943 AliFatal("Digits not found !") ;
944
945 //Take the digits list
946
85c25c2e 947 // Delete old if unzero
948 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
949 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
950 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
951 // Fill TRU and SM matrix
952 fAmpTrus = new TClonesArray("TMatrixD",nTRU);
953 fAmpTrus->SetName("AmpTrus");
954 fTimeRtrus = new TClonesArray("TMatrixD",nTRU);
955 fTimeRtrus->SetName("TimeRtrus");
956 fAmpSMods = new TClonesArray("TMatrixD",nSuperModules);
957 fAmpSMods->SetName("AmpSMods");
0964c2e9 958
85c25c2e 959 FillTRU(fDigitsList, fAmpTrus, fAmpSMods, fTimeRtrus);
960
0293e426 961 // Jet stuff - only one case, no freedom here
85c25c2e 962 if(fGeom->GetNEtaSubOfTRU() == 6) {
963 if(fAmpJetMatrix) {delete fAmpJetMatrix; fAmpJetMatrix=0;}
964 if(fJetMatrixE) {delete fJetMatrixE; fJetMatrixE=0;}
965
966 fAmpJetMatrix = new TMatrixD(17,12); // 17-phi(row), 12-eta(col)
967 fJetMatrixE = new TH2F("fJetMatrixE"," E of max patch in (#phi,#eta)",
968 17, 80.*TMath::DegToRad(), (180.+20.*2/3.)*TMath::DegToRad(), 12, -0.7, 0.7);
969 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row++) {
970 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
971 (*fAmpJetMatrix)(row,col) = 0.;
972 }
973 }
974 FillJetMatrixFromSMs(fAmpSMods, fAmpJetMatrix, fGeom);
975 }
976 if(!CheckConsistentOfMatrixes()) assert(0);
59264fa6 977
85c25c2e 978 // Do Tower Sliding and select Trigger
979 // Initialize varible that will contain maximum amplitudes and
980 // its corresponding tower position in eta and phi, and time.
981 TMatrixD ampmax2(4,nTRU) ; // 0-max amp, 1-irow, 2-icol, 3-timeR
9946f2fe 982 TMatrixD ampmaxn(4,nTRU) ;
0964c2e9 983
33d0b833 984 for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
0b2ec9f7 985 //Do 2x2 and nxn sums, select maximums.
85c25c2e 986
987 MakeSlidingTowers(fAmpTrus, fTimeRtrus, iSM, ampmax2, ampmaxn);
0964c2e9 988
85c25c2e 989 // Set the trigger
990 if(fIsolateInSuperModule) // here some discripency between tru and SM
991 SetTriggers(fAmpSMods,iSM,ampmax2,ampmaxn) ;
0964c2e9 992 if(!fIsolateInSuperModule)
85c25c2e 993 SetTriggers(fAmpTrus,iSM,ampmax2,ampmaxn) ;
59264fa6 994 }
0964c2e9 995
85c25c2e 996 // Do patch sliding and select Jet Trigger
997 // 0-max amp-meanFromVZERO(if), 1-irow, 2-icol, 3-timeR,
998 // 4-max amp , 5-meanFromVZERO (Nov 25, 2007)
999 // fAmpJetMax(6,1)
1000 MakeSlidingPatch((*fAmpJetMatrix), fNJetPatchPhi, fAmpJetMax); // no timing information here
1001
0964c2e9 1002 //Print();
85c25c2e 1003 // fDigitsList = 0;
1004}
1005
0293e426 1006//____________________________________________________________________________
581ee158 1007void AliEMCALTrigger::GetTriggerInfo(TArrayF &triggerPosition, TArrayF &triggerAmplitudes) const
85c25c2e 1008{
1009 // Template - should be defined; Nov 5, 2007
1010 triggerPosition[0] = 0.;
1011 triggerAmplitudes[0] = 0.;
1012}
1013
0293e426 1014//____________________________________________________________________________
85c25c2e 1015void AliEMCALTrigger::FillJetMatrixFromSMs(TClonesArray *ampmatrixsmod, TMatrixD* jetMat, AliEMCALGeometry *g)
1016{
1017 // Nov 5, 2007
1018 // Fill matrix for jet trigger from SM matrixes of modules
1019 //
1020 static int keyPrint = 0;
1021
1022 if(ampmatrixsmod==0 || jetMat==0 || g==0) return;
1023 Double_t amp = 0.0, ampSum=0.0;
0964c2e9 1024
85c25c2e 1025 Int_t nEtaModSum = g->GetNZ() / g->GetNEtaSubOfTRU(); // should be 4
1026 Int_t nPhiModSum = g->GetNPhi() / g->GetNTRUPhi(); // should be 4
1027
0293e426 1028 if(keyPrint) AliDebug(2,Form("%s",Form(" AliEMCALTrigger::FillJetMatrixFromSMs | nEtaModSum %i : nPhiModSum %i \n", nEtaModSum, nPhiModSum)));
85c25c2e 1029 Int_t jrow=0, jcol=0; // indexes of jet matrix
1030 Int_t nEtaSM=0, nPhiSM=0;
1031 for(Int_t iSM=0; iSM<ampmatrixsmod->GetEntries(); iSM++) {
1032 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSM));
a51e676d 1033
1034 if(!ampsmods) return;
1035
85c25c2e 1036 Int_t nrow = ampsmods->GetNrows();
1037 Int_t ncol = ampsmods->GetNcols();
1038 //printf("%s",Form(" ######## SM %i : nrow %i : ncol %i ##### \n", iSM, nrow, ncol));
1039 for(Int_t row=0; row<nrow; row++) {
1040 for(Int_t col=0; col<ncol; col++) {
1041 amp = (*ampsmods)(row,col);
1042 nPhiSM = iSM / 2;
1043 nEtaSM = iSM % 2;
1044 if (amp>0.0) {
0293e426 1045 if(keyPrint) AliDebug(2,Form("%s",Form(" ** nPhiSm %i : nEtaSM %i : row %2.2i : col %2.2i -> ", nPhiSM, nEtaSM, row, col)));
85c25c2e 1046 if(nEtaSM == 0) { // positive Z
1047 jrow = 3*nPhiSM + row/nPhiModSum;
1048 jcol = 6 + col / nEtaModSum;
1049 } else { // negative Z
1050 if(iSM<=9) jrow = 3*nPhiSM + 2 - row/nPhiModSum;
1051 else jrow = 3*nPhiSM + 1 - row/nPhiModSum; // half size
1052 jcol = 5 - col / nEtaModSum;
1053 }
0293e426 1054 if(keyPrint) AliDebug(2,Form("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat) \n", jrow, jcol, amp)));
85c25c2e 1055
1056 (*jetMat)(jrow,jcol) += amp;
1057 ampSum += amp; // For controling
1058 } else if(amp<0.0) {
8b518399 1059 AliDebug(1,Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat: amp<0) \n", jrow, jcol, amp));
85c25c2e 1060 assert(0);
1061 }
1062 }
1063 }
1064 } // cycle on SM
29b7e56e 1065 if(ampSum <= 0.0) AliDebug(1,Form("ampSum %f (<=0.0) ", ampSum));
85c25c2e 1066}
1067
0293e426 1068//____________________________________________________________________________
85c25c2e 1069void AliEMCALTrigger::MakeSlidingPatch(const TMatrixD &jm, const Int_t nPatchSize, TMatrixD &ampJetMax)
1070{
1071 // Sliding patch : nPatchSize x nPatchSize (OVERLAP)
1072 static int keyPrint = 0;
0293e426 1073 if(keyPrint) AliDebug(2,Form(" AliEMCALTrigger::MakeSlidingPatch() was started \n"));
85c25c2e 1074 Double_t ampCur = 0.0, e=0.0;
1075 ampJetMax(0,0) = 0.0;
1076 ampJetMax(3,0) = 0.0; // unused now
1077 ampJetMax(4,0) = ampJetMax(5,0) = 0.0;
1078 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row ++) {
1079 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
1080 ampCur = 0.;
1081 // check on patch size
1082 if( (row+nPatchSize-1) < fAmpJetMatrix->GetNrows() && (col+nPatchSize-1) < fAmpJetMatrix->GetNcols()){
1083 for(Int_t i = 0 ; i < nPatchSize ; i++) {
1084 for(Int_t j = 0 ; j < nPatchSize ; j++) {
1085 ampCur += jm(row+i, col+j);
1086 }
1087 } // end cycle on patch
1088 if(ampCur > ampJetMax(0,0)){
1089 ampJetMax(0,0) = ampCur;
1090 ampJetMax(1,0) = row;
1091 ampJetMax(2,0) = col;
1092 }
1093 } // check on patch size
1094 }
1095 }
0293e426 1096 if(keyPrint) AliDebug(2,Form(" ampJetMax %i row %2i->%2i col %2i->%2i \n", Int_t(ampJetMax(0,0)), Int_t(ampJetMax(1,0)), Int_t(ampJetMax(1,0))+nPatchSize-1, Int_t(ampJetMax(2,0)), Int_t(ampJetMax(2,0))+nPatchSize-1));
85c25c2e 1097
1098 Double_t eCorrJetMatrix=0.0;
1099 if(fVZER0Mult > 0.0) {
1100 // Correct patch energy (adc) and jet patch matrix energy
1101 Double_t meanAmpBG = GetMeanEmcalPatchEnergy(Int_t(fVZER0Mult), nPatchSize)/0.0153;
1102 ampJetMax(4,0) = ampJetMax(0,0);
1103 ampJetMax(5,0) = meanAmpBG;
1104
1105 Double_t eCorr = ampJetMax(0,0) - meanAmpBG;
0293e426 1106 AliDebug(2,Form(" ampJetMax(0,0) %f meanAmpBG %f eCorr %f : ampJetMax(4,0) %f \n",
1107 ampJetMax(0,0), meanAmpBG, eCorr, ampJetMax(5,0)));
85c25c2e 1108 ampJetMax(0,0) = eCorr;
1109 // --
1110 eCorrJetMatrix = GetMeanEmcalEnergy(Int_t(fVZER0Mult)) / 208.;
1111 }
1112 // Fill patch energy matrix
1113 for(int row=Int_t(ampJetMax(1,0)); row<Int_t(ampJetMax(1,0))+nPatchSize; row++) {
1114 for(int col=Int_t(ampJetMax(2,0)); col<Int_t(ampJetMax(2,0))+nPatchSize; col++) {
1115 e = Double_t(jm(row,col)*0.0153); // 0.0153 - hard coded now
1116 if(eCorrJetMatrix > 0.0) { // BG subtraction case
1117 e -= eCorrJetMatrix;
1118 fJetMatrixE->SetBinContent(row+1, col+1, e);
1119 } else if(e > 0.0) {
1120 fJetMatrixE->SetBinContent(row+1, col+1, e);
1121 }
1122 }
1123 }
1124 // PrintJetMatrix();
1125 // Set the jet trigger(s), multiple threshold now, Nov 19,2007
1126 for(Int_t i=0; i<fNJetThreshold; i++ ) {
1127 if(ampJetMax(0,0) >= fL1JetThreshold[i]) {
7e42cd80 1128 SetInput(GetNameOfJetTrigger(i));
85c25c2e 1129 }
1130 }
1131}
1132
0293e426 1133//____________________________________________________________________________
85c25c2e 1134Double_t AliEMCALTrigger::GetEmcalSumAmp() const
1135{
1136 // Return sum of amplidutes from EMCal
1137 // Used calibration coefficeint for transition to energy
1138 return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
1139}
1140
0293e426 1141//____________________________________________________________________________
85c25c2e 1142void AliEMCALTrigger::PrintJetMatrix() const
1143{
1144 // fAmpJetMatrix : (17,12); // 17-phi(row), 12-eta(col)
1145 if(fAmpJetMatrix == 0) return;
1146
0293e426 1147 AliInfo(Form("\n #### jetMatrix : (%i,%i) ##### \n ",
1148 fAmpJetMatrix->GetNrows(), fAmpJetMatrix->GetNcols()));
85c25c2e 1149 PrintMatrix(*fAmpJetMatrix);
1150}
1151
0293e426 1152//____________________________________________________________________________
85c25c2e 1153void AliEMCALTrigger::PrintAmpTruMatrix(Int_t ind) const
1154{
581ee158 1155 // Print matrix with TRU patches
85c25c2e 1156 TMatrixD * tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1157 if(tru == 0) return;
0293e426 1158 AliInfo(Form("\n #### Amp TRU matrix(%i) : (%i,%i) ##### \n ",
1159 ind, tru->GetNrows(), tru->GetNcols()));
85c25c2e 1160 PrintMatrix(*tru);
1161}
1162
0293e426 1163//____________________________________________________________________________
85c25c2e 1164void AliEMCALTrigger::PrintAmpSmMatrix(Int_t ind) const
1165{
581ee158 1166 // Print matrix with SM amplitudes
1167 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(ind));
85c25c2e 1168 if(sm == 0) return;
0293e426 1169 AliInfo(Form("\n #### Amp SM matrix(%i) : (%i,%i) ##### \n ",
1170 ind, sm->GetNrows(), sm->GetNcols()));
85c25c2e 1171 PrintMatrix(*sm);
1172}
1173
0293e426 1174//____________________________________________________________________________
85c25c2e 1175void AliEMCALTrigger::PrintMatrix(const TMatrixD &mat) const
1176{
581ee158 1177 //Print matrix object
0293e426 1178 for(Int_t col=0; col<mat.GetNcols(); col++) AliInfo(Form(" %3i ", col));
1179 AliInfo(Form("\n -- \n"));
85c25c2e 1180 for(Int_t row=0; row<mat.GetNrows(); row++) {
0293e426 1181 AliInfo(Form(" row:%2i ", row));
85c25c2e 1182 for(Int_t col=0; col<mat.GetNcols(); col++) {
0293e426 1183 AliInfo(Form(" %4i", (Int_t)mat(row,col)));
85c25c2e 1184 }
0293e426 1185 AliInfo("\n");
85c25c2e 1186 }
1187}
1188
0293e426 1189//____________________________________________________________________________
85c25c2e 1190Bool_t AliEMCALTrigger::CheckConsistentOfMatrixes(const Int_t pri)
1191{
581ee158 1192 // Check consitency of matrices
85c25c2e 1193 Double_t sumSM = 0.0, smCur=0.0;
1194 Double_t sumTru=0.0, sumTruInSM = 0.0, truSum=0.0;
1195 // Bool_t key = kTRUE;
1196
1197 for(Int_t i=0; i<fAmpSMods->GetEntries(); i++) {
1198 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(i));
1199 if(sm) {
1200 smCur = sm->Sum();
1201 sumSM += smCur;
1202
1203 sumTruInSM = 0.0;
1204 for(Int_t itru=0; itru<3; itru++) { // Cycle on tru inside SM
1205 Int_t ind = 3*i + itru;
1206 TMatrixD *tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1207 if(tru) {
1208 truSum = tru->Sum();
1209 sumTruInSM += truSum;
1210 }
1211 }
1212 sumTru += sumTruInSM;
1213
1214 if(sumTruInSM != smCur) {
8b518399 1215 AliDebug(1,Form(" sm %i : smCur %f -> sumTruInSM %f \n", i, smCur, sumTruInSM));
85c25c2e 1216 return kFALSE;
1217 }
1218 }
1219 }
1220 Double_t sumJetMat = fAmpJetMatrix->Sum();
581ee158 1221 if(pri || TMath::Abs(sumSM-sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001)
8b518399 1222 AliDebug(1,Form(" sumSM %f : sumTru %f : sumJetMat %f \n", sumSM, sumTru, sumJetMat));
581ee158 1223 if(TMath::Abs(sumSM - sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001) return kFALSE;
85c25c2e 1224 else return kTRUE;
1225}
1226
0293e426 1227//____________________________________________________________________________
85c25c2e 1228void AliEMCALTrigger::Browse(TBrowser* b)
1229{
581ee158 1230 //Browse.
85c25c2e 1231 if(&fInputs) b->Add(&fInputs);
1232 if(fAmpTrus) b->Add(fAmpTrus);
1233 if(fTimeRtrus) b->Add(fTimeRtrus);
1234 if(fAmpSMods) b->Add(fAmpSMods);
1235 if(fAmpJetMatrix) b->Add(fAmpJetMatrix);
1236 if(fJetMatrixE) b->Add(fJetMatrixE);
1237 // if(c) b->Add(c);
f0377b23 1238}