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