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