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