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