]> 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
correct z shift for C side
[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"
f0377b23 63
64ClassImp(AliEMCALTrigger)
65
85c25c2e 66TString AliEMCALTrigger::fgNameOfJetTriggers("EMCALJetTriggerL1");
67
f0377b23 68//______________________________________________________________________
69AliEMCALTrigger::AliEMCALTrigger()
9946f2fe 70 : AliTriggerDetector(), fGeom(0),
85c25c2e 71 f2x2MaxAmp(-1), f2x2ModulePhi(-1), f2x2ModuleEta(-1),
18a21c7c 72 f2x2SM(0),
85c25c2e 73 fnxnMaxAmp(-1), fnxnModulePhi(-1), fnxnModuleEta(-1),
0b2ec9f7 74 fnxnSM(0),
0964c2e9 75 fADCValuesHighnxn(0),fADCValuesLownxn(0),
76 fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
77 fDigitsList(0),
85c25c2e 78 fL0Threshold(100),fL1GammaLowPtThreshold(200),
79 fL1GammaMediumPtThreshold(500), fL1GammaHighPtThreshold(1000),
0964c2e9 80 fPatchSize(1), fIsolPatchSize(1),
81 f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1),
82 f2x2AmpOutOfPatchThres(100000), fnxnAmpOutOfPatchThres(100000),
83 fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),
85c25c2e 84 fSimulation(kTRUE), fIsolateInSuperModule(kTRUE), fTimeKey(kFALSE),
85 fAmpTrus(0),fTimeRtrus(0),fAmpSMods(0),
86 fTriggerPosition(6), fTriggerAmplitudes(4),
87 fNJetPatchPhi(3), fNJetPatchEta(3), fNJetThreshold(3), fL1JetThreshold(0), fJetMaxAmp(0),
88 fAmpJetMatrix(0), fJetMatrixE(0), fAmpJetMax(6,1), fVZER0Mult(0.)
f0377b23 89{
59264fa6 90 //ctor
0964c2e9 91 fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins];
92 fADCValuesLownxn = 0x0; //new Int_t[fTimeBins];
93 fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins];
94 fADCValuesLow2x2 = 0x0; //new Int_t[fTimeBins];
59264fa6 95
59264fa6 96 SetName("EMCAL");
85c25c2e 97 // Define jet threshold - can not change from outside now
98 // fNJetThreshold = 7; // For MB Pythia suppression
99 fNJetThreshold = 10; // Hijing
100 fL1JetThreshold = new Double_t[fNJetThreshold];
101 if(fNJetThreshold == 7) {
102 fL1JetThreshold[0] = 5./0.0153;
103 fL1JetThreshold[1] = 8./0.0153;
104 fL1JetThreshold[2] = 10./0.0153;
105 fL1JetThreshold[3] = 12./0.0153;
106 fL1JetThreshold[4] = 13./0.0153;
107 fL1JetThreshold[5] = 14./0.0153;
108 fL1JetThreshold[6] = 15./0.0153;
109 } else if(fNJetThreshold == 10) {
110 Double_t thGev[10]={5.,8.,10., 12., 13.,14.,15., 17., 20., 25.};
111 for(Int_t i=0; i<fNJetThreshold; i++) fL1JetThreshold[i] = thGev[i]/0.0153;
112 } else {
113 fL1JetThreshold[0] = 5./0.0153;
114 fL1JetThreshold[1] = 10./0.0153;
115 fL1JetThreshold[2] = 15./0.0153;
116 fL1JetThreshold[3] = 20./0.0153;
117 fL1JetThreshold[4] = 25./0.0153;
118 }
119 //
59264fa6 120 CreateInputs();
85c25c2e 121
122 fInputs.SetName("TriggersInputs");
59264fa6 123 //Print("") ;
f0377b23 124}
125
126
127
128//____________________________________________________________________________
129AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig)
18a21c7c 130 : AliTriggerDetector(trig),
9946f2fe 131 fGeom(trig.fGeom),
18a21c7c 132 f2x2MaxAmp(trig.f2x2MaxAmp),
85c25c2e 133 f2x2ModulePhi(trig.f2x2ModulePhi),
134 f2x2ModuleEta(trig.f2x2ModuleEta),
18a21c7c 135 f2x2SM(trig.f2x2SM),
0b2ec9f7 136 fnxnMaxAmp(trig.fnxnMaxAmp),
85c25c2e 137 fnxnModulePhi(trig.fnxnModulePhi),
138 fnxnModuleEta(trig.fnxnModuleEta),
0b2ec9f7 139 fnxnSM(trig.fnxnSM),
140 fADCValuesHighnxn(trig.fADCValuesHighnxn),
141 fADCValuesLownxn(trig.fADCValuesLownxn),
18a21c7c 142 fADCValuesHigh2x2(trig.fADCValuesHigh2x2),
143 fADCValuesLow2x2(trig.fADCValuesLow2x2),
144 fDigitsList(trig.fDigitsList),
145 fL0Threshold(trig.fL0Threshold),
85c25c2e 146 fL1GammaLowPtThreshold(trig.fL1GammaLowPtThreshold),
147 fL1GammaMediumPtThreshold(trig.fL1GammaMediumPtThreshold),
148 fL1GammaHighPtThreshold(trig.fL1GammaHighPtThreshold),
0964c2e9 149 fPatchSize(trig.fPatchSize),
150 fIsolPatchSize(trig.fIsolPatchSize),
151 f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch),
152 fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch),
153 f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres),
154 fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres),
155 fIs2x2Isol(trig.fIs2x2Isol),
156 fIsnxnIsol(trig.fIsnxnIsol),
157 fSimulation(trig.fSimulation),
85c25c2e 158 fIsolateInSuperModule(trig.fIsolateInSuperModule),
159 fTimeKey(trig.fTimeKey),
160 fAmpTrus(trig.fAmpTrus),
161 fTimeRtrus(trig.fTimeRtrus),
162 fAmpSMods(trig.fAmpSMods),
163 fTriggerPosition(trig.fTriggerPosition),
164 fTriggerAmplitudes(trig.fTriggerAmplitudes),
165 fNJetPatchPhi(trig.fNJetPatchPhi),
166 fNJetPatchEta(trig.fNJetPatchEta),
167 fNJetThreshold(trig.fNJetThreshold),
168 fL1JetThreshold(trig.fL1JetThreshold),
169 fJetMaxAmp(trig.fJetMaxAmp),
170 fAmpJetMatrix(trig.fAmpJetMatrix),
171 fJetMatrixE(trig.fJetMatrixE),
172 fAmpJetMax(trig.fAmpJetMax),
173 fVZER0Mult(trig.fVZER0Mult)
f0377b23 174{
f0377b23 175 // cpy ctor
f0377b23 176}
177
0293e426 178//____________________________________________________________________________
c35bbfd4 179AliEMCALTrigger::~AliEMCALTrigger() {
581ee158 180
181 //dtor
182
85c25c2e 183 if(GetTimeKey()) {
184 delete [] fADCValuesHighnxn;
185 delete [] fADCValuesLownxn;
186 delete [] fADCValuesHigh2x2;
187 delete [] fADCValuesLow2x2;
188 }
189 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
190 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
191 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
192 if(fAmpJetMatrix) delete fAmpJetMatrix;
193 if(fJetMatrixE) delete fJetMatrixE;
194 if(fL1JetThreshold) delete [] fL1JetThreshold;
c35bbfd4 195}
196
f0377b23 197//----------------------------------------------------------------------
198void AliEMCALTrigger::CreateInputs()
199{
200 // inputs
201
202 // Do not create inputs again!!
203 if( fInputs.GetEntriesFast() > 0 ) return;
85c25c2e 204
205 // Second parameter should be detector name = "EMCAL"
206 TString det("EMCAL"); // Apr 29, 2008
207 fInputs.AddLast( new AliTriggerInput( det+"_L0", det, 0x02) );
208 fInputs.AddLast( new AliTriggerInput( det+"_GammaHPt_L1", det, 0x04 ) );
209 fInputs.AddLast( new AliTriggerInput( det+"_GammaMPt_L1", det, 0x08 ) );
210 fInputs.AddLast( new AliTriggerInput( det+"_GammaLPt_L1", det, 0x016 ) );
ad534ff6 211 fInputs.AddLast( new AliTriggerInput( det+"_JetHPt_L1", det, 0x032 ) );
2994ff3b 212 fInputs.AddLast( new AliTriggerInput( det+"_JetMPt_L1", det, 0x048 ) );
213 fInputs.AddLast( new AliTriggerInput( det+"_JetLPt_L1", det, 0x064 ) );
85c25c2e 214
215 if(fNJetThreshold<=0) return;
216 // Jet Trigger(s)
217 UInt_t level = 0x032;
218 for(Int_t i=0; i<fNJetThreshold; i++ ) {
7e42cd80 219 TString name(GetNameOfJetTrigger(i));
85c25c2e 220 TString title("EMCAL Jet triger L1 :"); // unused now
221 // 0.0153 - hard coded now
222 title += Form("Th %i(%5.1f GeV) :", (Int_t)fL1JetThreshold[i], fL1JetThreshold[i] * 0.0153);
223 title += Form("patch %ix%i~(%3.2f(#phi)x%3.2f(#eta)) ",
224 fNJetPatchPhi, fNJetPatchEta, 0.11*(fNJetPatchPhi), 0.11*(fNJetPatchEta));
225 fInputs.AddLast( new AliTriggerInput(name, det, UChar_t(level)) );
226 level *= 2;
227 }
f0377b23 228
229}
230
231//____________________________________________________________________________
0964c2e9 232Bool_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) {
233
85c25c2e 234 // Nov 8, 2007
235 // EMCAL RTU size is 4modules(phi) x 24modules (eta)
236 // So maximum size of patch is 4modules x 4modules (EMCAL L0 trigger).
237 // Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch,
238 // inside isolation patch . iPatchType = 0 means calculation for 2x2 patch,
239 // iPatchType = 1 means calculation for nxn patch.
240 // In the next table there is an example of the different options of patch size and isolation patch size:
241 // Patch Size (fPatchSize)
242 // 0 1
243 // fIsolPatchSize 0 2x2 (not overlap) 4x4 (overlapped)
244 // 1 4x4 8x8
0964c2e9 245
246 Bool_t b = kFALSE;
0964c2e9 247
85c25c2e 248 // Get matrix of TRU or Module with maximum amplitude patch.
249 Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36.
0964c2e9 250 TMatrixD * ampmatrix = 0x0;
251 Int_t colborder = 0;
252 Int_t rowborder = 0;
85c25c2e 253 static int keyPrint = 0;
0293e426 254 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated : iSM %i mtru %i itru %i maxphi %i maxeta %i \n", iSM, mtru, itru, maxphi, maxeta));
0964c2e9 255
85c25c2e 256 if(fIsolateInSuperModule){ // ?
0964c2e9 257 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
85c25c2e 258 rowborder = fGeom->GetNPhi();
259 colborder = fGeom->GetNZ();
0964c2e9 260 AliDebug(2,"Isolate trigger in Module");
85c25c2e 261 } else{
0964c2e9 262 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
85c25c2e 263 rowborder = fGeom->GetNModulesInTRUPhi();
264 colborder = fGeom->GetNModulesInTRUEta();
0964c2e9 265 AliDebug(2,"Isolate trigger in TRU");
266 }
85c25c2e 267 if(iSM>9) rowborder /= 2; // half size in phi
0964c2e9 268
85c25c2e 269 //Define patch modules - what is this ??
270 Int_t isolmodules = fIsolPatchSize*(1+iPatchType);
271 Int_t ipatchmodules = 2*(1+fPatchSize*iPatchType);
272 Int_t minrow = maxphi - isolmodules;
273 Int_t mincol = maxeta - isolmodules;
274 Int_t maxrow = maxphi + isolmodules + ipatchmodules;
275 Int_t maxcol = maxeta + isolmodules + ipatchmodules;
276
277 minrow = minrow<0?0 :minrow;
278 mincol = mincol<0?0 :mincol;
279
280 maxrow = maxrow>rowborder?rowborder :maxrow;
281 maxcol = maxcol>colborder?colborder :maxcol;
9946f2fe 282
85c25c2e 283 //printf("%s\n",Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
284 //printf("%s\n",Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
285 // AliDebug(2,Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
286 //AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
0964c2e9 287
0964c2e9 288 //Add amplitudes in all isolation patch
85c25c2e 289 Float_t amp = 0.;
0964c2e9 290 for(Int_t irow = minrow ; irow < maxrow; irow ++)
291 for(Int_t icol = mincol ; icol < maxcol ; icol ++)
292 amp += (*ampmatrix)(irow,icol);
293
294 AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
295
296 if(amp < maxamp){
85c25c2e 297 // AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
298 // ampmatrix->Print();
0964c2e9 299 return kFALSE;
85c25c2e 300 } else {
0964c2e9 301 amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
85c25c2e 302 }
0964c2e9 303
304 AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
305
306 //Fill isolation amplitude data member and say if patch is isolated.
307 if(iPatchType == 0){ //2x2 case
308 f2x2AmpOutOfPatch = amp;
85c25c2e 309 if(amp < f2x2AmpOutOfPatchThres) b=kTRUE;
310 } else if(iPatchType == 1){ //nxn case
0964c2e9 311 fnxnAmpOutOfPatch = amp;
85c25c2e 312 if(amp < fnxnAmpOutOfPatchThres) b=kTRUE;
0964c2e9 313 }
314
0293e426 315 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated - OUT \n"));
85c25c2e 316
0964c2e9 317 return b;
318
319}
320
85c25c2e 321/*
0964c2e9 322//____________________________________________________________________________
c35bbfd4 323void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD &ampmax2, TMatrixD &ampmaxn){
f0377b23 324
85c25c2e 325 //Sums energy of all possible 2x2 (L0) and nxn (L1) modules per each TRU.
326 //Fast signal in the experiment is given by 2x2 modules,
327 //for this reason we loop inside the TRU modules by 2.
33d0b833 328
59264fa6 329 //Declare and initialize variables
9946f2fe 330 Int_t nCellsPhi = fGeom->GetNCellsInTRUPhi();
33d0b833 331 if(isupermod > 9)
59264fa6 332 nCellsPhi = nCellsPhi / 2 ; //Half size SM. Not Final.
f0377b23 333 // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
9946f2fe 334 Int_t nCellsEta = fGeom->GetNCellsInTRUEta();
335 Int_t nTRU = fGeom->GetNTRU();
f0377b23 336 // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
0964c2e9 337 //Int_t nTRU = geom->GeNTRU();//3 TRU per super module
f0377b23 338
59264fa6 339 Float_t amp2 = 0 ;
0b2ec9f7 340 Float_t ampn = 0 ;
33d0b833 341 for(Int_t i = 0; i < 4; i++){
9946f2fe 342 for(Int_t j = 0; j < nTRU; j++){
c35bbfd4 343 ampmax2(i,j) = -1;
344 ampmaxn(i,j) = -1;
59264fa6 345 }
346 }
f0377b23 347
59264fa6 348 //Create matrix that will contain 2x2 amplitude sums
0b2ec9f7 349 //used to calculate the nxn sums
c35bbfd4 350 TMatrixD tru2x2(nCellsPhi/2,nCellsEta/2) ;
59264fa6 351 for(Int_t i = 0; i < nCellsPhi/2; i++)
352 for(Int_t j = 0; j < nCellsEta/2; j++)
c35bbfd4 353 tru2x2(i,j) = -1;
59264fa6 354
355 //Loop over all TRUS in a supermodule
9946f2fe 356 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
59264fa6 357 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
358 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
9946f2fe 359 Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
33d0b833 360
59264fa6 361 //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
362 for(Int_t irow = 0 ; irow < nCellsPhi; irow += 2){
363 for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
364 amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
365 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
33d0b833 366
0964c2e9 367 //Fill matrix with added 2x2 cells for use in nxn sums
c35bbfd4 368 tru2x2(irow/2,icol/2) = amp2 ;
59264fa6 369 //Select 2x2 maximum sums to select L0
c35bbfd4 370 if(amp2 > ampmax2(0,mtru)){
371 ampmax2(0,mtru) = amp2 ;
372 ampmax2(1,mtru) = irow;
373 ampmax2(2,mtru) = icol;
59264fa6 374 }
375 }
376 }
377
378 //Find most recent time in the selected 2x2 cell
c35bbfd4 379 ampmax2(3,mtru) = 1 ;
380 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
381 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
59264fa6 382 for(Int_t i = 0; i<2; i++){
383 for(Int_t j = 0; j<2; j++){
384 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
c35bbfd4 385 if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) )
386 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j);
59264fa6 387 }
388 }
389 }
0b2ec9f7 390
391 //Sliding nxn, add nxn amplitudes (OVERLAP)
392 if(fPatchSize > 0){
393 for(Int_t irow = 0 ; irow < nCellsPhi/2; irow++){
394 for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
395 ampn = 0;
396 if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
397 for(Int_t i = 0 ; i <= fPatchSize ; i++)
398 for(Int_t j = 0 ; j <= fPatchSize ; j++)
c35bbfd4 399 ampn += tru2x2(irow+i,icol+j);
0b2ec9f7 400 //Select nxn maximum sums to select L1
c35bbfd4 401 if(ampn > ampmaxn(0,mtru)){
402 ampmaxn(0,mtru) = ampn ;
403 ampmaxn(1,mtru) = irow*2;
404 ampmaxn(2,mtru) = icol*2;
0b2ec9f7 405 }
59264fa6 406 }
407 }
408 }
0b2ec9f7 409
410 //Find most recent time in selected nxn cell
c35bbfd4 411 ampmaxn(3,mtru) = 1 ;
412 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
413 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
0b2ec9f7 414 for(Int_t i = 0; i<4*fPatchSize; i++){
415 for(Int_t j = 0; j<4*fPatchSize; j++){
9946f2fe 416 if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta){//Avoid exit the TRU
0b2ec9f7 417 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
c35bbfd4 418 if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) )
419 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j);
0b2ec9f7 420 }
421 }
59264fa6 422 }
423 }
424 }
0b2ec9f7 425 else {
c35bbfd4 426 ampmaxn(0,mtru) = ampmax2(0,mtru);
427 ampmaxn(1,mtru) = ampmax2(1,mtru);
428 ampmaxn(2,mtru) = ampmax2(2,mtru);
429 ampmaxn(3,mtru) = ampmax2(3,mtru);
0b2ec9f7 430 }
f0377b23 431 }
f0377b23 432}
85c25c2e 433*/
434//____________________________________________________________________________
435void AliEMCALTrigger::MakeSlidingTowers(const TClonesArray * amptrus, const TClonesArray * timeRtrus,
436const Int_t isupermod,TMatrixD &ampmax2, TMatrixD &ampmaxn){
437
438 // Output from module (2x2 cells from one module)
439 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi(); // now 4 modules (3 div in phi)
440 if(isupermod > 9)
441 nModulesPhi = nModulesPhi / 2 ; // Half size SM. Not Final.
442 //
443 Int_t nModulesEta = fGeom->GetNModulesInTRUEta(); // now 24 modules (no division in eta)
444 Int_t nTRU = fGeom->GetNTRU();
445 static int keyPrint = 0;
0293e426 446 if(keyPrint) AliDebug(2,Form("MakeSlidingTowers : nTRU %i nModulesPhi %i nModulesEta %i ",
447 nTRU, nModulesPhi, nModulesEta ));
85c25c2e 448
449 Float_t amp2 = 0 ;
450 Float_t ampn = 0 ;
451 for(Int_t i = 0; i < 4; i++){
452 for(Int_t j = 0; j < nTRU; j++){
453 ampmax2(i,j) = ampmaxn(i,j) = -1;
454 }
455 }
456
457 // Create matrix that will contain 2x2 amplitude sums
458 // used to calculate the nxn sums
459 TMatrixD tru2x2(nModulesPhi/2,nModulesEta/2);
460
461 // Loop over all TRUS in a supermodule
462 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
463 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
464 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
465 Int_t mtru = itru - isupermod*nTRU ; // Number of TRU in Supermodule !!
466
467 // Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
468 for(Int_t irow = 0 ; irow < nModulesPhi; irow +=2){
469 for(Int_t icol = 0 ; icol < nModulesEta ; icol +=2){
470 amp2 = (*amptru)(irow,icol) +(*amptru)(irow+1,icol)+
471 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
472
473 //Fill matrix with added 2x2 towers for use in nxn sums
474 tru2x2(irow/2,icol/2) = amp2 ;
475 //Select 2x2 maximum sums to select L0
476 if(amp2 > ampmax2(0,mtru)){
477 ampmax2(0,mtru) = amp2 ;
478 ampmax2(1,mtru) = irow;
479 ampmax2(2,mtru) = icol;
480 }
481 }
482 }
483
484 ampmax2(3,mtru) = 0.;
485 if(GetTimeKey()) {
486 // Find most recent time in the selected 2x2 towers
487 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
488 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
489 for(Int_t i = 0; i<2; i++){
490 for(Int_t j = 0; j<2; j++){
491 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
492 if((*timeRtru)(row2+i,col2+j) > ampmax2(3,mtru) )
493 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); // max time
494 }
495 }
496 }
497 }
498
499 //Sliding nxn, add nxn amplitudes (OVERLAP)
500 if(fPatchSize > 0){
501 for(Int_t irow = 0 ; irow < nModulesPhi/2; irow++){
502 for(Int_t icol = 0 ; icol < nModulesEta/2; icol++){
503 ampn = 0;
504 if( (irow+fPatchSize) < nModulesPhi/2 && (icol+fPatchSize) < nModulesEta/2){ //Avoid exit the TRU
505 for(Int_t i = 0 ; i <= fPatchSize ; i++)
506 for(Int_t j = 0 ; j <= fPatchSize ; j++)
507 ampn += tru2x2(irow+i,icol+j);
508 //Select nxn maximum sums to select L1
509 if(ampn > ampmaxn(0,mtru)){
510 ampmaxn(0,mtru) = ampn ;
511 ampmaxn(1,mtru) = irow;
512 ampmaxn(2,mtru) = icol;
513 }
514 }
515 }
516 }
517
518 ampmaxn(3,mtru) = 0.; // Was 1 , I don't know why
519 if(GetTimeKey()) {
520 //Find most recent time in selected nxn cell
521 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
522 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
523 for(Int_t i = 0; i<4*fPatchSize; i++){
524 for(Int_t j = 0; j<4*fPatchSize; j++){
525 if( (rown+i) < nModulesPhi && (coln+j) < nModulesEta){//Avoid exit the TRU
526 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
527 if((*timeRtru)(rown+i,coln+j) > ampmaxn(3,mtru) )
528 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); // max time
529 }
530 }
531 }
532 }
533 }
534 } else { // copy 2x2 to nxn
535 ampmaxn(0,mtru) = ampmax2(0,mtru);
536 ampmaxn(1,mtru) = ampmax2(1,mtru);
537 ampmaxn(2,mtru) = ampmax2(2,mtru);
538 ampmaxn(3,mtru) = ampmax2(3,mtru);
539 }
540 }
0293e426 541 if(keyPrint) AliDebug(2,Form(" : MakeSlidingTowers -OUt \n"));
85c25c2e 542}
f0377b23 543
544//____________________________________________________________________________
545void AliEMCALTrigger::Print(const Option_t * opt) const
546{
547
548 //Prints main parameters
549
550 if(! opt)
551 return;
552 AliTriggerInput* in = 0x0 ;
0293e426 553 AliInfo(Form(" fSimulation %i (input option) : #digits %i\n", fSimulation, fDigitsList->GetEntries()));
554 AliInfo(Form(" fTimeKey %i \n ", fTimeKey));
555
556 AliInfo(Form("\t Maximum Amplitude after Sliding Cell, \n")) ;
557 AliInfo(Form("\t -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
558 f2x2MaxAmp,f2x2SM)) ;
559 AliInfo(Form("\t -2x2 from row %d to row %d and from column %d to column %d\n", f2x2ModulePhi, f2x2ModulePhi+2, f2x2ModuleEta, f2x2ModuleEta+2));
560 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 561 if(fPatchSize > 0){
0293e426 562 AliInfo(Form("\t Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1)));
563 AliInfo(Form("\t -nxn cells sum (overlapped) : %10.2f, in Super Module %d\n", fnxnMaxAmp,fnxnSM));
564 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)) ;
565 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 566 }
0964c2e9 567
0293e426 568 AliInfo(Form("\t Isolate in SuperModule? %d\n", fIsolateInSuperModule)) ;
569 AliInfo(Form("\t Threshold for LO %10.2f\n", fL0Threshold));
0964c2e9 570
59264fa6 571 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
f0377b23 572 if(in->GetValue())
0293e426 573 AliInfo(Form("\t *** EMCAL LO is set ***\n"));
f0377b23 574
0293e426 575 AliInfo(Form("\t Gamma Low Pt Threshold for L1 %10.2f\n", fL1GammaLowPtThreshold));
85c25c2e 576 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_GammaLPt_L1" );
f0377b23 577 if(in->GetValue())
0293e426 578 AliInfo(Form("\t *** EMCAL Gamma Low Pt for L1 is set ***\n"));
f0377b23 579
0293e426 580 AliInfo(Form("\t Gamma Medium Pt Threshold for L1 %10.2f\n", fL1GammaMediumPtThreshold));
85c25c2e 581 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaMPt_L1" );
f0377b23 582 if(in->GetValue())
0293e426 583 AliInfo(Form("\t *** EMCAL Gamma Medium Pt for L1 is set ***\n"));
f0377b23 584
0293e426 585 AliInfo(Form("\t Gamma High Pt Threshold for L1 %10.2f\n", fL1GammaHighPtThreshold));
85c25c2e 586 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaHPt_L1" );
f0377b23 587 if(in->GetValue())
0293e426 588 AliInfo(Form("\t *** 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)) ;
829ba234 815 amp = Float_t(dig->GetAmplitude()); // Energy of the digit (arbitrary units)
85c25c2e 816 id = dig->GetId() ; // Id label of the cell
817 timeR = dig->GetTimeR() ; // Earliest time of the digit
8b518399 818 if(amp<=0.0) AliDebug(1,Form(" id %i amp %f \n", id, amp));
85c25c2e 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)
0293e426 823 AliError(Form("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}
0293e426 869
59264fa6 870//____________________________________________________________________________
871void AliEMCALTrigger::Trigger()
872{
873 //Main Method to select triggers.
581ee158 874 TH1::AddDirectory(0);
875
33c3c91a 876 AliRunLoader *runLoader = AliRunLoader::Instance();
85c25c2e 877 AliEMCALLoader *emcalLoader = 0;
878 if(runLoader) {
879 emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
880 }
0964c2e9 881
59264fa6 882 //Load EMCAL Geometry
85c25c2e 883 if (runLoader && runLoader->GetAliRun() && runLoader->GetAliRun()->GetDetector("EMCAL"))
9946f2fe 884 fGeom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
0964c2e9 885
384c0bba 886 if (fGeom == 0)
887 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
888
9946f2fe 889 if (fGeom==0)
59264fa6 890 AliFatal("Did not get geometry from EMCALLoader");
0964c2e9 891
59264fa6 892 //Define parameters
9946f2fe 893 Int_t nSuperModules = fGeom->GetNumberOfSuperModules() ; //12 SM in EMCAL
85c25c2e 894 Int_t nTRU = fGeom->GetNTRU(); // 3 TRU per super module
59264fa6 895
896 //Intialize data members each time the trigger is called in event loop
85c25c2e 897 f2x2MaxAmp = -1; f2x2ModulePhi = -1; f2x2ModuleEta = -1;
898 fnxnMaxAmp = -1; fnxnModulePhi = -1; fnxnModuleEta = -1;
59264fa6 899
85c25c2e 900 // Take the digits list if simulation
901 if(fSimulation && runLoader && emcalLoader){ // works than run seperate macros
c787fb51 902 runLoader->LoadDigits("EMCAL");
59264fa6 903 fDigitsList = emcalLoader->Digits() ;
85c25c2e 904 runLoader->LoadSDigits("EMCAL");
59264fa6 905 }
85c25c2e 906 // Digits list should be set by method SetDigitsList(TClonesArray * digits)
59264fa6 907 if(!fDigitsList)
908 AliFatal("Digits not found !") ;
909
910 //Take the digits list
911
85c25c2e 912 // Delete old if unzero
913 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
914 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
915 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
916 // Fill TRU and SM matrix
917 fAmpTrus = new TClonesArray("TMatrixD",nTRU);
918 fAmpTrus->SetName("AmpTrus");
919 fTimeRtrus = new TClonesArray("TMatrixD",nTRU);
920 fTimeRtrus->SetName("TimeRtrus");
921 fAmpSMods = new TClonesArray("TMatrixD",nSuperModules);
922 fAmpSMods->SetName("AmpSMods");
0964c2e9 923
85c25c2e 924 FillTRU(fDigitsList, fAmpTrus, fAmpSMods, fTimeRtrus);
925
0293e426 926 // Jet stuff - only one case, no freedom here
85c25c2e 927 if(fGeom->GetNEtaSubOfTRU() == 6) {
928 if(fAmpJetMatrix) {delete fAmpJetMatrix; fAmpJetMatrix=0;}
929 if(fJetMatrixE) {delete fJetMatrixE; fJetMatrixE=0;}
930
931 fAmpJetMatrix = new TMatrixD(17,12); // 17-phi(row), 12-eta(col)
932 fJetMatrixE = new TH2F("fJetMatrixE"," E of max patch in (#phi,#eta)",
933 17, 80.*TMath::DegToRad(), (180.+20.*2/3.)*TMath::DegToRad(), 12, -0.7, 0.7);
934 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row++) {
935 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
936 (*fAmpJetMatrix)(row,col) = 0.;
937 }
938 }
939 FillJetMatrixFromSMs(fAmpSMods, fAmpJetMatrix, fGeom);
940 }
941 if(!CheckConsistentOfMatrixes()) assert(0);
59264fa6 942
85c25c2e 943 // Do Tower Sliding and select Trigger
944 // Initialize varible that will contain maximum amplitudes and
945 // its corresponding tower position in eta and phi, and time.
946 TMatrixD ampmax2(4,nTRU) ; // 0-max amp, 1-irow, 2-icol, 3-timeR
9946f2fe 947 TMatrixD ampmaxn(4,nTRU) ;
0964c2e9 948
33d0b833 949 for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
0b2ec9f7 950 //Do 2x2 and nxn sums, select maximums.
85c25c2e 951
952 MakeSlidingTowers(fAmpTrus, fTimeRtrus, iSM, ampmax2, ampmaxn);
0964c2e9 953
85c25c2e 954 // Set the trigger
955 if(fIsolateInSuperModule) // here some discripency between tru and SM
956 SetTriggers(fAmpSMods,iSM,ampmax2,ampmaxn) ;
0964c2e9 957 if(!fIsolateInSuperModule)
85c25c2e 958 SetTriggers(fAmpTrus,iSM,ampmax2,ampmaxn) ;
59264fa6 959 }
0964c2e9 960
85c25c2e 961 // Do patch sliding and select Jet Trigger
962 // 0-max amp-meanFromVZERO(if), 1-irow, 2-icol, 3-timeR,
963 // 4-max amp , 5-meanFromVZERO (Nov 25, 2007)
964 // fAmpJetMax(6,1)
965 MakeSlidingPatch((*fAmpJetMatrix), fNJetPatchPhi, fAmpJetMax); // no timing information here
966
0964c2e9 967 //Print();
85c25c2e 968 // fDigitsList = 0;
969}
970
0293e426 971//____________________________________________________________________________
581ee158 972void AliEMCALTrigger::GetTriggerInfo(TArrayF &triggerPosition, TArrayF &triggerAmplitudes) const
85c25c2e 973{
974 // Template - should be defined; Nov 5, 2007
975 triggerPosition[0] = 0.;
976 triggerAmplitudes[0] = 0.;
977}
978
0293e426 979//____________________________________________________________________________
85c25c2e 980void AliEMCALTrigger::FillJetMatrixFromSMs(TClonesArray *ampmatrixsmod, TMatrixD* jetMat, AliEMCALGeometry *g)
981{
982 // Nov 5, 2007
983 // Fill matrix for jet trigger from SM matrixes of modules
984 //
985 static int keyPrint = 0;
986
987 if(ampmatrixsmod==0 || jetMat==0 || g==0) return;
988 Double_t amp = 0.0, ampSum=0.0;
0964c2e9 989
85c25c2e 990 Int_t nEtaModSum = g->GetNZ() / g->GetNEtaSubOfTRU(); // should be 4
991 Int_t nPhiModSum = g->GetNPhi() / g->GetNTRUPhi(); // should be 4
992
0293e426 993 if(keyPrint) AliDebug(2,Form("%s",Form(" AliEMCALTrigger::FillJetMatrixFromSMs | nEtaModSum %i : nPhiModSum %i \n", nEtaModSum, nPhiModSum)));
85c25c2e 994 Int_t jrow=0, jcol=0; // indexes of jet matrix
995 Int_t nEtaSM=0, nPhiSM=0;
996 for(Int_t iSM=0; iSM<ampmatrixsmod->GetEntries(); iSM++) {
997 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSM));
998 Int_t nrow = ampsmods->GetNrows();
999 Int_t ncol = ampsmods->GetNcols();
1000 //printf("%s",Form(" ######## SM %i : nrow %i : ncol %i ##### \n", iSM, nrow, ncol));
1001 for(Int_t row=0; row<nrow; row++) {
1002 for(Int_t col=0; col<ncol; col++) {
1003 amp = (*ampsmods)(row,col);
1004 nPhiSM = iSM / 2;
1005 nEtaSM = iSM % 2;
1006 if (amp>0.0) {
0293e426 1007 if(keyPrint) AliDebug(2,Form("%s",Form(" ** nPhiSm %i : nEtaSM %i : row %2.2i : col %2.2i -> ", nPhiSM, nEtaSM, row, col)));
85c25c2e 1008 if(nEtaSM == 0) { // positive Z
1009 jrow = 3*nPhiSM + row/nPhiModSum;
1010 jcol = 6 + col / nEtaModSum;
1011 } else { // negative Z
1012 if(iSM<=9) jrow = 3*nPhiSM + 2 - row/nPhiModSum;
1013 else jrow = 3*nPhiSM + 1 - row/nPhiModSum; // half size
1014 jcol = 5 - col / nEtaModSum;
1015 }
0293e426 1016 if(keyPrint) AliDebug(2,Form("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat) \n", jrow, jcol, amp)));
85c25c2e 1017
1018 (*jetMat)(jrow,jcol) += amp;
1019 ampSum += amp; // For controling
1020 } else if(amp<0.0) {
8b518399 1021 AliDebug(1,Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat: amp<0) \n", jrow, jcol, amp));
85c25c2e 1022 assert(0);
1023 }
1024 }
1025 }
1026 } // cycle on SM
9310749e 1027 if(ampSum <= 0.0) AliDebug(1,Form("FillJetMatrixFromSMs","ampSum %f (<=0.0) ", ampSum));
85c25c2e 1028}
1029
0293e426 1030//____________________________________________________________________________
85c25c2e 1031void AliEMCALTrigger::MakeSlidingPatch(const TMatrixD &jm, const Int_t nPatchSize, TMatrixD &ampJetMax)
1032{
1033 // Sliding patch : nPatchSize x nPatchSize (OVERLAP)
1034 static int keyPrint = 0;
0293e426 1035 if(keyPrint) AliDebug(2,Form(" AliEMCALTrigger::MakeSlidingPatch() was started \n"));
85c25c2e 1036 Double_t ampCur = 0.0, e=0.0;
1037 ampJetMax(0,0) = 0.0;
1038 ampJetMax(3,0) = 0.0; // unused now
1039 ampJetMax(4,0) = ampJetMax(5,0) = 0.0;
1040 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row ++) {
1041 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
1042 ampCur = 0.;
1043 // check on patch size
1044 if( (row+nPatchSize-1) < fAmpJetMatrix->GetNrows() && (col+nPatchSize-1) < fAmpJetMatrix->GetNcols()){
1045 for(Int_t i = 0 ; i < nPatchSize ; i++) {
1046 for(Int_t j = 0 ; j < nPatchSize ; j++) {
1047 ampCur += jm(row+i, col+j);
1048 }
1049 } // end cycle on patch
1050 if(ampCur > ampJetMax(0,0)){
1051 ampJetMax(0,0) = ampCur;
1052 ampJetMax(1,0) = row;
1053 ampJetMax(2,0) = col;
1054 }
1055 } // check on patch size
1056 }
1057 }
0293e426 1058 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 1059
1060 Double_t eCorrJetMatrix=0.0;
1061 if(fVZER0Mult > 0.0) {
1062 // Correct patch energy (adc) and jet patch matrix energy
1063 Double_t meanAmpBG = GetMeanEmcalPatchEnergy(Int_t(fVZER0Mult), nPatchSize)/0.0153;
1064 ampJetMax(4,0) = ampJetMax(0,0);
1065 ampJetMax(5,0) = meanAmpBG;
1066
1067 Double_t eCorr = ampJetMax(0,0) - meanAmpBG;
0293e426 1068 AliDebug(2,Form(" ampJetMax(0,0) %f meanAmpBG %f eCorr %f : ampJetMax(4,0) %f \n",
1069 ampJetMax(0,0), meanAmpBG, eCorr, ampJetMax(5,0)));
85c25c2e 1070 ampJetMax(0,0) = eCorr;
1071 // --
1072 eCorrJetMatrix = GetMeanEmcalEnergy(Int_t(fVZER0Mult)) / 208.;
1073 }
1074 // Fill patch energy matrix
1075 for(int row=Int_t(ampJetMax(1,0)); row<Int_t(ampJetMax(1,0))+nPatchSize; row++) {
1076 for(int col=Int_t(ampJetMax(2,0)); col<Int_t(ampJetMax(2,0))+nPatchSize; col++) {
1077 e = Double_t(jm(row,col)*0.0153); // 0.0153 - hard coded now
1078 if(eCorrJetMatrix > 0.0) { // BG subtraction case
1079 e -= eCorrJetMatrix;
1080 fJetMatrixE->SetBinContent(row+1, col+1, e);
1081 } else if(e > 0.0) {
1082 fJetMatrixE->SetBinContent(row+1, col+1, e);
1083 }
1084 }
1085 }
1086 // PrintJetMatrix();
1087 // Set the jet trigger(s), multiple threshold now, Nov 19,2007
1088 for(Int_t i=0; i<fNJetThreshold; i++ ) {
1089 if(ampJetMax(0,0) >= fL1JetThreshold[i]) {
7e42cd80 1090 SetInput(GetNameOfJetTrigger(i));
85c25c2e 1091 }
1092 }
1093}
1094
0293e426 1095//____________________________________________________________________________
85c25c2e 1096Double_t AliEMCALTrigger::GetEmcalSumAmp() const
1097{
1098 // Return sum of amplidutes from EMCal
1099 // Used calibration coefficeint for transition to energy
1100 return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
1101}
1102
0293e426 1103//____________________________________________________________________________
85c25c2e 1104void AliEMCALTrigger::PrintJetMatrix() const
1105{
1106 // fAmpJetMatrix : (17,12); // 17-phi(row), 12-eta(col)
1107 if(fAmpJetMatrix == 0) return;
1108
0293e426 1109 AliInfo(Form("\n #### jetMatrix : (%i,%i) ##### \n ",
1110 fAmpJetMatrix->GetNrows(), fAmpJetMatrix->GetNcols()));
85c25c2e 1111 PrintMatrix(*fAmpJetMatrix);
1112}
1113
0293e426 1114//____________________________________________________________________________
85c25c2e 1115void AliEMCALTrigger::PrintAmpTruMatrix(Int_t ind) const
1116{
581ee158 1117 // Print matrix with TRU patches
85c25c2e 1118 TMatrixD * tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1119 if(tru == 0) return;
0293e426 1120 AliInfo(Form("\n #### Amp TRU matrix(%i) : (%i,%i) ##### \n ",
1121 ind, tru->GetNrows(), tru->GetNcols()));
85c25c2e 1122 PrintMatrix(*tru);
1123}
1124
0293e426 1125//____________________________________________________________________________
85c25c2e 1126void AliEMCALTrigger::PrintAmpSmMatrix(Int_t ind) const
1127{
581ee158 1128 // Print matrix with SM amplitudes
1129 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(ind));
85c25c2e 1130 if(sm == 0) return;
0293e426 1131 AliInfo(Form("\n #### Amp SM matrix(%i) : (%i,%i) ##### \n ",
1132 ind, sm->GetNrows(), sm->GetNcols()));
85c25c2e 1133 PrintMatrix(*sm);
1134}
1135
0293e426 1136//____________________________________________________________________________
85c25c2e 1137void AliEMCALTrigger::PrintMatrix(const TMatrixD &mat) const
1138{
581ee158 1139 //Print matrix object
0293e426 1140 for(Int_t col=0; col<mat.GetNcols(); col++) AliInfo(Form(" %3i ", col));
1141 AliInfo(Form("\n -- \n"));
85c25c2e 1142 for(Int_t row=0; row<mat.GetNrows(); row++) {
0293e426 1143 AliInfo(Form(" row:%2i ", row));
85c25c2e 1144 for(Int_t col=0; col<mat.GetNcols(); col++) {
0293e426 1145 AliInfo(Form(" %4i", (Int_t)mat(row,col)));
85c25c2e 1146 }
0293e426 1147 AliInfo("\n");
85c25c2e 1148 }
1149}
1150
0293e426 1151//____________________________________________________________________________
85c25c2e 1152Bool_t AliEMCALTrigger::CheckConsistentOfMatrixes(const Int_t pri)
1153{
581ee158 1154 // Check consitency of matrices
85c25c2e 1155 Double_t sumSM = 0.0, smCur=0.0;
1156 Double_t sumTru=0.0, sumTruInSM = 0.0, truSum=0.0;
1157 // Bool_t key = kTRUE;
1158
1159 for(Int_t i=0; i<fAmpSMods->GetEntries(); i++) {
1160 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(i));
1161 if(sm) {
1162 smCur = sm->Sum();
1163 sumSM += smCur;
1164
1165 sumTruInSM = 0.0;
1166 for(Int_t itru=0; itru<3; itru++) { // Cycle on tru inside SM
1167 Int_t ind = 3*i + itru;
1168 TMatrixD *tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1169 if(tru) {
1170 truSum = tru->Sum();
1171 sumTruInSM += truSum;
1172 }
1173 }
1174 sumTru += sumTruInSM;
1175
1176 if(sumTruInSM != smCur) {
8b518399 1177 AliDebug(1,Form(" sm %i : smCur %f -> sumTruInSM %f \n", i, smCur, sumTruInSM));
85c25c2e 1178 return kFALSE;
1179 }
1180 }
1181 }
1182 Double_t sumJetMat = fAmpJetMatrix->Sum();
581ee158 1183 if(pri || TMath::Abs(sumSM-sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001)
8b518399 1184 AliDebug(1,Form(" sumSM %f : sumTru %f : sumJetMat %f \n", sumSM, sumTru, sumJetMat));
581ee158 1185 if(TMath::Abs(sumSM - sumTru)>0.0001 || TMath::Abs(sumSM-sumJetMat) > 0.0001) return kFALSE;
85c25c2e 1186 else return kTRUE;
1187}
1188
0293e426 1189//____________________________________________________________________________
85c25c2e 1190void AliEMCALTrigger::Browse(TBrowser* b)
1191{
581ee158 1192 //Browse.
85c25c2e 1193 if(&fInputs) b->Add(&fInputs);
1194 if(fAmpTrus) b->Add(fAmpTrus);
1195 if(fTimeRtrus) b->Add(fTimeRtrus);
1196 if(fAmpSMods) b->Add(fAmpSMods);
1197 if(fAmpJetMatrix) b->Add(fAmpJetMatrix);
1198 if(fJetMatrixE) b->Add(fJetMatrixE);
1199 // if(c) b->Add(c);
f0377b23 1200}