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