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