1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //_________________________________________________________________________
19 // Class for trigger analysis.
20 // Digits are grouped in TRU's (Trigger Units). A TRU consists of 384
21 // modules ordered fNTRUPhi x fNTRUEta. The algorithm searches all possible 2x2
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.
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.
31 // //Inside the event loop
32 // AliEMCALTrigger *tr = new AliEMCALTrigger();//Init Trigger
33 // tr->SetL0Threshold(100); //Arbitrary threshold values
34 // tr->SetL1GammaLowPtThreshold(1000);
35 // tr->SetL1GammaMediumPtThreshold(10000);
36 // tr->SetL1GammaHighPtThreshold(20000);
38 // tr->Trigger(); //Execute Trigger
39 // tr->Print(""); //Print results
41 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN)
42 //////////////////////////////////////////////////////////////////////////////
46 // --- ROOT system ---
52 // --- ALIROOT system ---
54 #include "AliRunLoader.h"
55 #include "AliTriggerInput.h"
57 #include "AliEMCALLoader.h"
58 #include "AliEMCALDigit.h"
59 #include "AliEMCALTrigger.h"
60 #include "AliEMCALGeometry.h"
61 #include "AliEMCALRawUtils.h"
64 ClassImp(AliEMCALTrigger)
66 TString AliEMCALTrigger::fgNameOfJetTriggers("EMCALJetTriggerL1");
68 //______________________________________________________________________
69 AliEMCALTrigger::AliEMCALTrigger()
70 : AliTriggerDetector(), fGeom(0),
71 f2x2MaxAmp(-1), f2x2ModulePhi(-1), f2x2ModuleEta(-1),
73 fnxnMaxAmp(-1), fnxnModulePhi(-1), fnxnModuleEta(-1),
75 fADCValuesHighnxn(0),fADCValuesLownxn(0),
76 fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
78 fL0Threshold(100),fL1GammaLowPtThreshold(200),
79 fL1GammaMediumPtThreshold(500), fL1GammaHighPtThreshold(1000),
80 fPatchSize(1), fIsolPatchSize(1),
81 f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1),
82 f2x2AmpOutOfPatchThres(100000), fnxnAmpOutOfPatchThres(100000),
83 fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),
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.)
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];
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;
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;
122 fInputs.SetName("TriggersInputs");
128 //____________________________________________________________________________
129 AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig)
130 : AliTriggerDetector(trig),
132 f2x2MaxAmp(trig.f2x2MaxAmp),
133 f2x2ModulePhi(trig.f2x2ModulePhi),
134 f2x2ModuleEta(trig.f2x2ModuleEta),
136 fnxnMaxAmp(trig.fnxnMaxAmp),
137 fnxnModulePhi(trig.fnxnModulePhi),
138 fnxnModuleEta(trig.fnxnModuleEta),
140 fADCValuesHighnxn(trig.fADCValuesHighnxn),
141 fADCValuesLownxn(trig.fADCValuesLownxn),
142 fADCValuesHigh2x2(trig.fADCValuesHigh2x2),
143 fADCValuesLow2x2(trig.fADCValuesLow2x2),
144 fDigitsList(trig.fDigitsList),
145 fL0Threshold(trig.fL0Threshold),
146 fL1GammaLowPtThreshold(trig.fL1GammaLowPtThreshold),
147 fL1GammaMediumPtThreshold(trig.fL1GammaMediumPtThreshold),
148 fL1GammaHighPtThreshold(trig.fL1GammaHighPtThreshold),
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),
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)
178 //____________________________________________________________________________
179 AliEMCALTrigger::~AliEMCALTrigger() {
181 delete [] fADCValuesHighnxn;
182 delete [] fADCValuesLownxn;
183 delete [] fADCValuesHigh2x2;
184 delete [] fADCValuesLow2x2;
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;
194 //----------------------------------------------------------------------
195 void AliEMCALTrigger::CreateInputs()
199 // Do not create inputs again!!
200 if( fInputs.GetEntriesFast() > 0 ) return;
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 fInputs.AddLast( new AliTriggerInput( det+"_JetHPt_L1", det, 0x032 ) );
209 fInputs.AddLast( new AliTriggerInput( det+"_JetMPt_L1", det, 0x048 ) );
210 fInputs.AddLast( new AliTriggerInput( det+"_JetLPt_L1", det, 0x064 ) );
212 if(fNJetThreshold<=0) return;
214 UInt_t level = 0x032;
215 for(Int_t i=0; i<fNJetThreshold; i++ ) {
216 TString name(GetNameOfJetTrigger(i));
217 TString title("EMCAL Jet triger L1 :"); // unused now
218 // 0.0153 - hard coded now
219 title += Form("Th %i(%5.1f GeV) :", (Int_t)fL1JetThreshold[i], fL1JetThreshold[i] * 0.0153);
220 title += Form("patch %ix%i~(%3.2f(#phi)x%3.2f(#eta)) ",
221 fNJetPatchPhi, fNJetPatchEta, 0.11*(fNJetPatchPhi), 0.11*(fNJetPatchEta));
222 fInputs.AddLast( new AliTriggerInput(name, det, UChar_t(level)) );
228 //____________________________________________________________________________
229 Bool_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) {
232 // EMCAL RTU size is 4modules(phi) x 24modules (eta)
233 // So maximum size of patch is 4modules x 4modules (EMCAL L0 trigger).
234 // Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch,
235 // inside isolation patch . iPatchType = 0 means calculation for 2x2 patch,
236 // iPatchType = 1 means calculation for nxn patch.
237 // In the next table there is an example of the different options of patch size and isolation patch size:
238 // Patch Size (fPatchSize)
240 // fIsolPatchSize 0 2x2 (not overlap) 4x4 (overlapped)
245 // Get matrix of TRU or Module with maximum amplitude patch.
246 Int_t itru = mtru + iSM * fGeom->GetNTRU(); //number of tru, min 0 max 3*12=36.
247 TMatrixD * ampmatrix = 0x0;
250 static int keyPrint = 0;
251 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated : iSM %i mtru %i itru %i maxphi %i maxeta %i \n", iSM, mtru, itru, maxphi, maxeta));
253 if(fIsolateInSuperModule){ // ?
254 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
255 rowborder = fGeom->GetNPhi();
256 colborder = fGeom->GetNZ();
257 AliDebug(2,"Isolate trigger in Module");
259 ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
260 rowborder = fGeom->GetNModulesInTRUPhi();
261 colborder = fGeom->GetNModulesInTRUEta();
262 AliDebug(2,"Isolate trigger in TRU");
264 if(iSM>9) rowborder /= 2; // half size in phi
266 //Define patch modules - what is this ??
267 Int_t isolmodules = fIsolPatchSize*(1+iPatchType);
268 Int_t ipatchmodules = 2*(1+fPatchSize*iPatchType);
269 Int_t minrow = maxphi - isolmodules;
270 Int_t mincol = maxeta - isolmodules;
271 Int_t maxrow = maxphi + isolmodules + ipatchmodules;
272 Int_t maxcol = maxeta + isolmodules + ipatchmodules;
274 minrow = minrow<0?0 :minrow;
275 mincol = mincol<0?0 :mincol;
277 maxrow = maxrow>rowborder?rowborder :maxrow;
278 maxcol = maxcol>colborder?colborder :maxcol;
280 //printf("%s\n",Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
281 //printf("%s\n",Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
282 // AliDebug(2,Form("Number of added Isol Modules %d, Patch Size %d",isolmodules, ipatchmodules));
283 //AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
285 //Add amplitudes in all isolation patch
287 for(Int_t irow = minrow ; irow < maxrow; irow ++)
288 for(Int_t icol = mincol ; icol < maxcol ; icol ++)
289 amp += (*ampmatrix)(irow,icol);
291 AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
294 // AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
295 // ampmatrix->Print();
298 amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
301 AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
303 //Fill isolation amplitude data member and say if patch is isolated.
304 if(iPatchType == 0){ //2x2 case
305 f2x2AmpOutOfPatch = amp;
306 if(amp < f2x2AmpOutOfPatchThres) b=kTRUE;
307 } else if(iPatchType == 1){ //nxn case
308 fnxnAmpOutOfPatch = amp;
309 if(amp < fnxnAmpOutOfPatchThres) b=kTRUE;
312 if(keyPrint) AliDebug(2,Form(" IsPatchIsolated - OUT \n"));
319 //____________________________________________________________________________
320 void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
322 //Sums energy of all possible 2x2 (L0) and nxn (L1) modules per each TRU.
323 //Fast signal in the experiment is given by 2x2 modules,
324 //for this reason we loop inside the TRU modules by 2.
326 //Declare and initialize variables
327 Int_t nCellsPhi = fGeom->GetNCellsInTRUPhi();
329 nCellsPhi = nCellsPhi / 2 ; //Half size SM. Not Final.
330 // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
331 Int_t nCellsEta = fGeom->GetNCellsInTRUEta();
332 Int_t nTRU = fGeom->GetNTRU();
333 // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
334 //Int_t nTRU = geom->GeNTRU();//3 TRU per super module
338 for(Int_t i = 0; i < 4; i++){
339 for(Int_t j = 0; j < nTRU; j++){
345 //Create matrix that will contain 2x2 amplitude sums
346 //used to calculate the nxn sums
347 TMatrixD tru2x2(nCellsPhi/2,nCellsEta/2) ;
348 for(Int_t i = 0; i < nCellsPhi/2; i++)
349 for(Int_t j = 0; j < nCellsEta/2; j++)
352 //Loop over all TRUS in a supermodule
353 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
354 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
355 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
356 Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
358 //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
359 for(Int_t irow = 0 ; irow < nCellsPhi; irow += 2){
360 for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
361 amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
362 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
364 //Fill matrix with added 2x2 cells for use in nxn sums
365 tru2x2(irow/2,icol/2) = amp2 ;
366 //Select 2x2 maximum sums to select L0
367 if(amp2 > ampmax2(0,mtru)){
368 ampmax2(0,mtru) = amp2 ;
369 ampmax2(1,mtru) = irow;
370 ampmax2(2,mtru) = icol;
375 //Find most recent time in the selected 2x2 cell
376 ampmax2(3,mtru) = 1 ;
377 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
378 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
379 for(Int_t i = 0; i<2; i++){
380 for(Int_t j = 0; j<2; j++){
381 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
382 if((*timeRtru)(row2+i,col2+j) < ampmax2(3,mtru) )
383 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j);
388 //Sliding nxn, add nxn amplitudes (OVERLAP)
390 for(Int_t irow = 0 ; irow < nCellsPhi/2; irow++){
391 for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
393 if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
394 for(Int_t i = 0 ; i <= fPatchSize ; i++)
395 for(Int_t j = 0 ; j <= fPatchSize ; j++)
396 ampn += tru2x2(irow+i,icol+j);
397 //Select nxn maximum sums to select L1
398 if(ampn > ampmaxn(0,mtru)){
399 ampmaxn(0,mtru) = ampn ;
400 ampmaxn(1,mtru) = irow*2;
401 ampmaxn(2,mtru) = icol*2;
407 //Find most recent time in selected nxn cell
408 ampmaxn(3,mtru) = 1 ;
409 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
410 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
411 for(Int_t i = 0; i<4*fPatchSize; i++){
412 for(Int_t j = 0; j<4*fPatchSize; j++){
413 if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta){//Avoid exit the TRU
414 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
415 if((*timeRtru)(rown+i,coln+j) < ampmaxn(3,mtru) )
416 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j);
423 ampmaxn(0,mtru) = ampmax2(0,mtru);
424 ampmaxn(1,mtru) = ampmax2(1,mtru);
425 ampmaxn(2,mtru) = ampmax2(2,mtru);
426 ampmaxn(3,mtru) = ampmax2(3,mtru);
431 //____________________________________________________________________________
432 void AliEMCALTrigger::MakeSlidingTowers(const TClonesArray * amptrus, const TClonesArray * timeRtrus,
433 const Int_t isupermod,TMatrixD &max2, TMatrixD &maxn){
435 // Output from module (2x2 cells from one module)
436 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi(); // now 4 modules (3 div in phi)
438 nModulesPhi = nModulesPhi / 2 ; // Half size SM. Not Final.
440 Int_t nModulesEta = fGeom->GetNModulesInTRUEta(); // now 24 modules (no division in eta)
441 Int_t nTRU = fGeom->GetNTRU();
442 static int keyPrint = 0;
443 if(keyPrint) AliDebug(2,Form("MakeSlidingTowers : nTRU %i nModulesPhi %i nModulesEta %i ",
444 nTRU, nModulesPhi, nModulesEta ));
448 for(Int_t i = 0; i < 4; i++){
449 for(Int_t j = 0; j < nTRU; j++){
450 ampmax2(i,j) = ampmaxn(i,j) = -1;
454 // Create matrix that will contain 2x2 amplitude sums
455 // used to calculate the nxn sums
456 TMatrixD tru2x2(nModulesPhi/2,nModulesEta/2);
458 // Loop over all TRUS in a supermodule
459 for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
460 TMatrixD * amptru = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
461 TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
462 Int_t mtru = itru - isupermod*nTRU ; // Number of TRU in Supermodule !!
464 // Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
465 for(Int_t irow = 0 ; irow < nModulesPhi; irow +=2){
466 for(Int_t icol = 0 ; icol < nModulesEta ; icol +=2){
467 amp2 = (*amptru)(irow,icol) +(*amptru)(irow+1,icol)+
468 (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
470 //Fill matrix with added 2x2 towers for use in nxn sums
471 tru2x2(irow/2,icol/2) = amp2 ;
472 //Select 2x2 maximum sums to select L0
473 if(amp2 > ampmax2(0,mtru)){
474 ampmax2(0,mtru) = amp2 ;
475 ampmax2(1,mtru) = irow;
476 ampmax2(2,mtru) = icol;
481 ampmax2(3,mtru) = 0.;
483 // Find most recent time in the selected 2x2 towers
484 Int_t row2 = static_cast <Int_t> (ampmax2(1,mtru));
485 Int_t col2 = static_cast <Int_t> (ampmax2(2,mtru));
486 for(Int_t i = 0; i<2; i++){
487 for(Int_t j = 0; j<2; j++){
488 if((*amptru)(row2+i,col2+j) > 0 && (*timeRtru)(row2+i,col2+j)> 0){
489 if((*timeRtru)(row2+i,col2+j) > ampmax2(3,mtru) )
490 ampmax2(3,mtru) = (*timeRtru)(row2+i,col2+j); // max time
496 //Sliding nxn, add nxn amplitudes (OVERLAP)
498 for(Int_t irow = 0 ; irow < nModulesPhi/2; irow++){
499 for(Int_t icol = 0 ; icol < nModulesEta/2; icol++){
501 if( (irow+fPatchSize) < nModulesPhi/2 && (icol+fPatchSize) < nModulesEta/2){ //Avoid exit the TRU
502 for(Int_t i = 0 ; i <= fPatchSize ; i++)
503 for(Int_t j = 0 ; j <= fPatchSize ; j++)
504 ampn += tru2x2(irow+i,icol+j);
505 //Select nxn maximum sums to select L1
506 if(ampn > ampmaxn(0,mtru)){
507 ampmaxn(0,mtru) = ampn ;
508 ampmaxn(1,mtru) = irow;
509 ampmaxn(2,mtru) = icol;
515 ampmaxn(3,mtru) = 0.; // Was 1 , I don't know why
517 //Find most recent time in selected nxn cell
518 Int_t rown = static_cast <Int_t> (ampmaxn(1,mtru));
519 Int_t coln = static_cast <Int_t> (ampmaxn(2,mtru));
520 for(Int_t i = 0; i<4*fPatchSize; i++){
521 for(Int_t j = 0; j<4*fPatchSize; j++){
522 if( (rown+i) < nModulesPhi && (coln+j) < nModulesEta){//Avoid exit the TRU
523 if((*amptru)(rown+i,coln+j) > 0 && (*timeRtru)(rown+i,coln+j)> 0){
524 if((*timeRtru)(rown+i,coln+j) > ampmaxn(3,mtru) )
525 ampmaxn(3,mtru) = (*timeRtru)(rown+i,coln+j); // max time
531 } else { // copy 2x2 to nxn
532 ampmaxn(0,mtru) = ampmax2(0,mtru);
533 ampmaxn(1,mtru) = ampmax2(1,mtru);
534 ampmaxn(2,mtru) = ampmax2(2,mtru);
535 ampmaxn(3,mtru) = ampmax2(3,mtru);
538 if(keyPrint) AliDebug(2,Form(" : MakeSlidingTowers -OUt \n"));
541 //____________________________________________________________________________
542 void AliEMCALTrigger::Print(const Option_t * opt) const
545 //Prints main parameters
549 AliTriggerInput* in = 0x0 ;
550 AliInfo(Form(" fSimulation %i (input option) : #digits %i\n", fSimulation, fDigitsList->GetEntries()));
551 AliInfo(Form(" fTimeKey %i \n ", fTimeKey));
553 AliInfo(Form("\t Maximum Amplitude after Sliding Cell, \n")) ;
554 AliInfo(Form("\t -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
555 f2x2MaxAmp,f2x2SM)) ;
556 AliInfo(Form("\t -2x2 from row %d to row %d and from column %d to column %d\n", f2x2ModulePhi, f2x2ModulePhi+2, f2x2ModuleEta, f2x2ModuleEta+2));
557 AliInfo(Form("\t -2x2 Isolation Patch %d x %d, Amplitude out of 2x2 patch is %f, threshold %f, Isolated? %d \n", 2*fIsolPatchSize+2, 2*fIsolPatchSize+2, f2x2AmpOutOfPatch, f2x2AmpOutOfPatchThres,static_cast<Int_t> (fIs2x2Isol)));
559 AliInfo(Form("\t Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1)));
560 AliInfo(Form("\t -nxn cells sum (overlapped) : %10.2f, in Super Module %d\n", fnxnMaxAmp,fnxnSM));
561 AliInfo(Form("\t -nxn from row %d to row %d and from column %d to column %d\n", fnxnModulePhi, fnxnModulePhi+4*fPatchSize, fnxnModuleEta, fnxnModuleEta+4*fPatchSize)) ;
562 AliInfo(Form("\t -nxn Isolation Patch %d x %d, Amplitude out of nxn patch is %f, threshold %f, Isolated? %d \n", 4*fIsolPatchSize+2*(fPatchSize+1),4*fIsolPatchSize+2*(fPatchSize+1) , fnxnAmpOutOfPatch, fnxnAmpOutOfPatchThres,static_cast<Int_t> (fIsnxnIsol) ));
565 AliInfo(Form("\t Isolate in SuperModule? %d\n", fIsolateInSuperModule)) ;
566 AliInfo(Form("\t Threshold for LO %10.2f\n", fL0Threshold));
568 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
570 AliInfo(Form("\t *** EMCAL LO is set ***\n"));
572 AliInfo(Form("\t Gamma Low Pt Threshold for L1 %10.2f\n", fL1GammaLowPtThreshold));
573 in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_GammaLPt_L1" );
575 AliInfo(Form("\t *** EMCAL Gamma Low Pt for L1 is set ***\n"));
577 AliInfo(Form("\t Gamma Medium Pt Threshold for L1 %10.2f\n", fL1GammaMediumPtThreshold));
578 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaMPt_L1" );
580 AliInfo(Form("\t *** EMCAL Gamma Medium Pt for L1 is set ***\n"));
582 AliInfo(Form("\t Gamma High Pt Threshold for L1 %10.2f\n", fL1GammaHighPtThreshold));
583 in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_GammaHPt_L1" );
585 AliInfo(Form("\t *** EMCAL Gamma High Pt for L1 is set ***\n")) ;
589 //____________________________________________________________________________
590 void AliEMCALTrigger::SetTriggers(const TClonesArray * ampmatrix,const Int_t iSM,
591 const TMatrixD &max2,
592 const TMatrixD &maxn)
594 //Checks the 2x2 and nxn maximum amplitude per each TRU and
595 //compares with the different L0 and L1 triggers thresholds
596 Float_t max2[] = {-1,-1,-1,-1} ;
597 Float_t maxn[] = {-1,-1,-1,-1} ;
601 Int_t nTRU = fGeom->GetNTRU();
603 //Find maximum summed amplitude of all the TRU
605 for(Int_t i = 0 ; i < nTRU ; i++){
606 if(max2[0] < ampmax2(0,i) ){
607 max2[0] = ampmax2(0,i) ; // 2x2 summed max amplitude
608 max2[1] = ampmax2(1,i) ; // corresponding phi position in TRU
609 max2[2] = ampmax2(2,i) ; // corresponding eta position in TRU
610 max2[3] = ampmax2(3,i) ; // corresponding most recent time
613 if(maxn[0] < ampmaxn(0,i) ){
614 maxn[0] = ampmaxn(0,i) ; // nxn summed max amplitude
615 maxn[1] = ampmaxn(1,i) ; // corresponding phi position in TRU
616 maxn[2] = ampmaxn(2,i) ; // corresponding eta position in TRU
617 maxn[3] = ampmaxn(3,i) ; // corresponding most recent time
622 //--------Set max amplitude if larger than in other Super Modules------------
623 Float_t maxtimeR2 = -1 ;
624 Float_t maxtimeRn = -1 ;
625 static AliEMCALRawUtils rawUtil;
626 Int_t nTimeBins = rawUtil.GetRawFormatTimeBins() ;
628 //Set max of 2x2 amplitudes and select L0 trigger
629 if(max2[0] > f2x2MaxAmp ){
630 // if(max2[0] > 5) printf(" L0 : iSM %i: max2[0] %5.0f : max2[3] %5.0f (maxtimeR2) \n",
631 // iSM, max2[0], max2[3]);
632 f2x2MaxAmp = max2[0] ;
634 maxtimeR2 = max2[3] ;
635 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtru2,
636 static_cast<Int_t>(max2[1]),
637 static_cast<Int_t>(max2[2]),
638 f2x2ModulePhi,f2x2ModuleEta);
641 if(fIsolateInSuperModule)
642 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, f2x2ModulePhi,f2x2ModuleEta) ;
644 fIs2x2Isol = IsPatchIsolated(0, ampmatrix, iSM, mtru2, f2x2MaxAmp, static_cast<Int_t>(max2[1]), static_cast<Int_t>(max2[2])) ;
647 //Transform digit amplitude in Raw Samples
648 if (fADCValuesLow2x2 == 0) {
649 fADCValuesLow2x2 = new Int_t[nTimeBins];
650 fADCValuesHigh2x2 = new Int_t[nTimeBins];
652 //printf(" maxtimeR2 %12.5e (1)\n", maxtimeR2);
653 rawUtil.RawSampledResponse(maxtimeR2 * AliEMCALRawUtils::GetRawFormatTimeBin(),
654 f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ;
656 // Set Trigger Inputs, compare ADC time bins until threshold is attained
658 for(Int_t i = 0 ; i < nTimeBins ; i++){
659 // printf(" fADCValuesHigh2x2[%i] %i : %i \n", i, fADCValuesHigh2x2[i], fADCValuesLow2x2[i]);
660 if(fADCValuesHigh2x2[i] >= fL0Threshold || fADCValuesLow2x2[i] >= fL0Threshold){
661 SetInput("EMCAL_L0") ;
666 // Nov 5 - no analysis of time information
667 if(f2x2MaxAmp >= fL0Threshold) { // should add the low amp too
668 SetInput("EMCAL_L0");
673 //------------Set max of nxn amplitudes and select L1 trigger---------
674 if(maxn[0] > fnxnMaxAmp ){
675 fnxnMaxAmp = maxn[0] ;
677 maxtimeRn = maxn[3] ;
678 fGeom->GetModulePhiEtaIndexInSModuleFromTRUIndex(mtrun,
679 static_cast<Int_t>(maxn[1]),
680 static_cast<Int_t>(maxn[2]),
681 fnxnModulePhi,fnxnModuleEta) ;
684 if(fIsolateInSuperModule)
685 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, fnxnModulePhi, fnxnModuleEta) ;
687 fIsnxnIsol = IsPatchIsolated(1, ampmatrix, iSM, mtrun, fnxnMaxAmp, static_cast<Int_t>(maxn[1]), static_cast<Int_t>(maxn[2])) ;
690 //Transform digit amplitude in Raw Samples
691 if (fADCValuesLownxn == 0) {
692 fADCValuesHighnxn = new Int_t[nTimeBins];
693 fADCValuesLownxn = new Int_t[nTimeBins];
695 rawUtil.RawSampledResponse(maxtimeRn * AliEMCALRawUtils::GetRawFormatTimeBin(),
696 fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
698 //Set Trigger Inputs, compare ADC time bins until threshold is attained
700 for(Int_t i = 0 ; i < nTimeBins ; i++){
701 if(fADCValuesHighnxn[i] >= fL1GammaLowPtThreshold || fADCValuesLownxn[i] >= fL1GammaLowPtThreshold){
702 SetInput("EMCAL_GammaLPt_L1") ;
708 for(Int_t i = 0 ; i < nTimeBins ; i++){
709 if(fADCValuesHighnxn[i] >= fL1GammaMediumPtThreshold || fADCValuesLownxn[i] >= fL1GammaMediumPtThreshold){
710 SetInput("EMCAL_GammaMPt_L1") ;
716 for(Int_t i = 0 ; i < nTimeBins ; i++){
717 if(fADCValuesHighnxn[i] >= fL1GammaHighPtThreshold || fADCValuesLownxn[i] >= fL1GammaHighPtThreshold){
718 SetInput("EMCAL_GammaHPt_L1") ;
723 // Nov 5 - no analysis of time information
724 if(fnxnMaxAmp >= fL1GammaLowPtThreshold) { // should add the low amp too
725 SetInput("EMCAL_GammaLPt_L1") ; //SetL1 Low
727 if(fnxnMaxAmp >= fL1GammaMediumPtThreshold) { // should add the low amp too
728 SetInput("EMCAL_GammaMPt_L1") ; //SetL1 Medium
730 if(fnxnMaxAmp >= fL1GammaHighPtThreshold) { // should add the low amp too
731 SetInput("EMCAL_GammaHPt_L1") ; //SetL1 High
737 //____________________________________________________________________________
738 void AliEMCALTrigger::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * ampmatrixsmod, TClonesArray * timeRmatrix) {
740 // Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule.
741 // Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of
742 // TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
743 // Last 2 modules are half size in Phi, I considered that the number of TRU
744 // is maintained for the last modules but decision not taken. If different,
745 // then this must be changed. Also fill a matrix with all amplitudes in supermodule for isolation studies.
747 // Initilize and declare variables
748 // List of TRU matrices initialized to 0.
749 // printf("<I> AliEMCALTrigger::FillTRU() started : # digits %i\n", digits->GetEntriesFast());
752 // One input per EMCAL module so size of matrix is reduced by 4 (2x2 division case)
754 Int_t nPhi = fGeom->GetNPhi();
755 Int_t nZ = fGeom->GetNZ();
756 Int_t nTRU = fGeom->GetNTRU();
757 // Int_t nTRUPhi = fGeom->GetNTRUPhi();
758 Int_t nModulesPhi = fGeom->GetNModulesInTRUPhi();
759 Int_t nModulesPhi2 = fGeom->GetNModulesInTRUPhi();
760 Int_t nModulesEta = fGeom->GetNModulesInTRUEta();
761 // printf("<I> AliEMCALTrigger::FillTRU() nTRU %i nTRUPhi %i : nModulesPhi %i nModulesEta %i \n",
762 // nTRU, nTRUPhi, nModulesPhi, nModulesEta);
773 // iphim, ietam - module indexes in SM
777 //List of TRU matrices initialized to 0.
778 Int_t nSup = fGeom->GetNumberOfSuperModules();
779 for(Int_t k = 0; k < nTRU*nSup; k++){
780 TMatrixD amptrus(nModulesPhi,nModulesEta) ;
781 TMatrixD timeRtrus(nModulesPhi,nModulesEta) ;
782 // Do we need to initialise? I think TMatrixD does it by itself...
783 for(Int_t i = 0; i < nModulesPhi; i++){
784 for(Int_t j = 0; j < nModulesEta; j++){
786 timeRtrus(i,j) = 0.0;
789 new((*ampmatrix)[k]) TMatrixD(amptrus) ;
790 new((*timeRmatrix)[k]) TMatrixD(timeRtrus) ;
793 // List of Modules matrices initialized to 0.
794 for(Int_t k = 0; k < nSup ; k++){
796 // if(nSup>9) mphi = nPhi/2; // the same size
797 TMatrixD ampsmods( mphi, nZ);
798 for(Int_t i = 0; i < mphi; i++){
799 for(Int_t j = 0; j < nZ; j++){
803 new((*ampmatrixsmod)[k]) TMatrixD(ampsmods) ;
806 AliEMCALDigit * dig ;
808 //Digits loop to fill TRU matrices with amplitudes.
809 for(Int_t idig = 0 ; idig < digits->GetEntriesFast() ; idig++){
811 dig = dynamic_cast<AliEMCALDigit *>(digits->At(idig)) ;
812 amp = Float_t(dig->GetAmp()); // Energy of the digit (arbitrary units)
813 id = dig->GetId() ; // Id label of the cell
814 timeR = dig->GetTimeR() ; // Earliest time of the digit
815 if(amp<=0.0) AliDebug(1,Form(" id %i amp %f \n", id, amp));
816 // printf(" FILLTRU : timeR %10.5e time %10.5e : amp %10.5e \n", timeR, dig->GetTime(), amp);
817 // Get eta and phi cell position in supermodule
818 Bool_t bCell = fGeom->GetCellIndex(id, iSupMod, nModule, nIphi, nIeta) ;
820 AliError(Form("FillTRU","%i Wrong cell id number %i ", idig, id)) ;
822 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
823 // iphim, ietam - module indexes in SM
824 fGeom->GetModuleIndexesFromCellIndexesInSModule(iSupMod,iphi,ieta, iphim, ietam, nModule);
826 //printf("iSupMod %i nModule %i iphi %i ieta %i iphim %i ietam %i \n",
827 //iSupMod,nModule, iphi, ieta, iphim, ietam);
829 // Check to which TRU in the supermodule belongs the cell.
830 // Supermodules are divided in a TRU matrix of dimension
831 // (fNTRUPhi,fNTRUEta).
832 // Each TRU is a cell matrix of dimension (nModulesPhi,nModulesEta)
834 // First calculate the row and column in the supermodule
835 // of the TRU to which the cell belongs.
836 Int_t row = iphim / nModulesPhi;
837 Int_t col = ietam / nModulesEta;
838 //Calculate label number of the TRU
839 Int_t itru = fGeom->GetAbsTRUNumberFromNumberInSm(row, col, iSupMod);
841 //Fill TRU matrix with cell values
842 TMatrixD * amptrus = dynamic_cast<TMatrixD *>(ampmatrix->At(itru)) ;
843 TMatrixD * timeRtrus = dynamic_cast<TMatrixD *>(timeRmatrix->At(itru)) ;
845 //Calculate row and column of the module inside the TRU with number itru
846 Int_t irow = iphim - row * nModulesPhi;
848 irow = iphim - row * nModulesPhi2; // size of matrix the same
849 Int_t icol = ietam - col * nModulesEta;
851 (*amptrus)(irow,icol) += amp ;
852 if((*timeRtrus)(irow,icol) <0.0 || (*timeRtrus)(irow,icol) <= timeR){ // ??
853 (*timeRtrus)(irow,icol) = timeR ;
855 //printf(" ieta %i iphi %i iSM %i || col %i row %i : itru %i -> amp %f\n",
856 // ieta, iphi, iSupMod, col, row, itru, amp);
857 //####################SUPERMODULE MATRIX ##################
858 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSupMod)) ;
859 (*ampsmods)(iphim,ietam) += amp ;
860 // printf(" id %i iphim %i ietam %i SM %i : irow %i icol %i itru %i : amp %6.0f\n",
861 //id, iphim, ietam, iSupMod, irow, icol, itru, amp);
864 //printf("<I> AliEMCALTrigger::FillTRU() is ended \n");
867 //____________________________________________________________________________
868 void AliEMCALTrigger::Trigger()
870 TH1::AddDirectory(0);
871 //Main Method to select triggers.
872 AliRunLoader *runLoader = AliRunLoader::Instance();
873 AliEMCALLoader *emcalLoader = 0;
875 emcalLoader = dynamic_cast<AliEMCALLoader*>(runLoader->GetDetectorLoader("EMCAL"));
878 //Load EMCAL Geometry
879 if (runLoader && runLoader->GetAliRun() && runLoader->GetAliRun()->GetDetector("EMCAL"))
880 fGeom = dynamic_cast<AliEMCAL*>(runLoader->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
883 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
886 AliFatal("Did not get geometry from EMCALLoader");
889 Int_t nSuperModules = fGeom->GetNumberOfSuperModules() ; //12 SM in EMCAL
890 Int_t nTRU = fGeom->GetNTRU(); // 3 TRU per super module
892 //Intialize data members each time the trigger is called in event loop
893 f2x2MaxAmp = -1; f2x2ModulePhi = -1; f2x2ModuleEta = -1;
894 fnxnMaxAmp = -1; fnxnModulePhi = -1; fnxnModuleEta = -1;
896 // Take the digits list if simulation
897 if(fSimulation && runLoader && emcalLoader){ // works than run seperate macros
898 runLoader->LoadDigits("EMCAL");
899 fDigitsList = emcalLoader->Digits() ;
900 runLoader->LoadSDigits("EMCAL");
902 // Digits list should be set by method SetDigitsList(TClonesArray * digits)
904 AliFatal("Digits not found !") ;
906 //Take the digits list
908 // Delete old if unzero
909 if(fAmpTrus) {fAmpTrus->Delete(); delete fAmpTrus;}
910 if(fTimeRtrus) {fTimeRtrus->Delete(); delete fTimeRtrus;}
911 if(fAmpSMods) {fAmpSMods->Delete(); delete fAmpSMods;}
912 // Fill TRU and SM matrix
913 fAmpTrus = new TClonesArray("TMatrixD",nTRU);
914 fAmpTrus->SetName("AmpTrus");
915 fTimeRtrus = new TClonesArray("TMatrixD",nTRU);
916 fTimeRtrus->SetName("TimeRtrus");
917 fAmpSMods = new TClonesArray("TMatrixD",nSuperModules);
918 fAmpSMods->SetName("AmpSMods");
920 FillTRU(fDigitsList, fAmpTrus, fAmpSMods, fTimeRtrus);
922 // Jet stuff - only one case, no freedom here
923 if(fGeom->GetNEtaSubOfTRU() == 6) {
924 if(fAmpJetMatrix) {delete fAmpJetMatrix; fAmpJetMatrix=0;}
925 if(fJetMatrixE) {delete fJetMatrixE; fJetMatrixE=0;}
927 fAmpJetMatrix = new TMatrixD(17,12); // 17-phi(row), 12-eta(col)
928 fJetMatrixE = new TH2F("fJetMatrixE"," E of max patch in (#phi,#eta)",
929 17, 80.*TMath::DegToRad(), (180.+20.*2/3.)*TMath::DegToRad(), 12, -0.7, 0.7);
930 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row++) {
931 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
932 (*fAmpJetMatrix)(row,col) = 0.;
935 FillJetMatrixFromSMs(fAmpSMods, fAmpJetMatrix, fGeom);
937 if(!CheckConsistentOfMatrixes()) assert(0);
939 // Do Tower Sliding and select Trigger
940 // Initialize varible that will contain maximum amplitudes and
941 // its corresponding tower position in eta and phi, and time.
942 TMatrixD ampmax2(4,nTRU) ; // 0-max amp, 1-irow, 2-icol, 3-timeR
943 TMatrixD ampmaxn(4,nTRU) ;
945 for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
946 //Do 2x2 and nxn sums, select maximums.
948 MakeSlidingTowers(fAmpTrus, fTimeRtrus, iSM, ampmax2, ampmaxn);
951 if(fIsolateInSuperModule) // here some discripency between tru and SM
952 SetTriggers(fAmpSMods,iSM,ampmax2,ampmaxn) ;
953 if(!fIsolateInSuperModule)
954 SetTriggers(fAmpTrus,iSM,ampmax2,ampmaxn) ;
957 // Do patch sliding and select Jet Trigger
958 // 0-max amp-meanFromVZERO(if), 1-irow, 2-icol, 3-timeR,
959 // 4-max amp , 5-meanFromVZERO (Nov 25, 2007)
961 MakeSlidingPatch((*fAmpJetMatrix), fNJetPatchPhi, fAmpJetMax); // no timing information here
967 //____________________________________________________________________________
968 void AliEMCALTrigger::GetTriggerInfo(TArrayF &triggerPosition, TArrayF &triggerAmplitudes)
970 // Template - should be defined; Nov 5, 2007
971 triggerPosition[0] = 0.;
972 triggerAmplitudes[0] = 0.;
975 //____________________________________________________________________________
976 void AliEMCALTrigger::FillJetMatrixFromSMs(TClonesArray *ampmatrixsmod, TMatrixD* jetMat, AliEMCALGeometry *g)
979 // Fill matrix for jet trigger from SM matrixes of modules
981 static int keyPrint = 0;
983 if(ampmatrixsmod==0 || jetMat==0 || g==0) return;
984 Double_t amp = 0.0, ampSum=0.0;
986 Int_t nEtaModSum = g->GetNZ() / g->GetNEtaSubOfTRU(); // should be 4
987 Int_t nPhiModSum = g->GetNPhi() / g->GetNTRUPhi(); // should be 4
989 if(keyPrint) AliDebug(2,Form("%s",Form(" AliEMCALTrigger::FillJetMatrixFromSMs | nEtaModSum %i : nPhiModSum %i \n", nEtaModSum, nPhiModSum)));
990 Int_t jrow=0, jcol=0; // indexes of jet matrix
991 Int_t nEtaSM=0, nPhiSM=0;
992 for(Int_t iSM=0; iSM<ampmatrixsmod->GetEntries(); iSM++) {
993 TMatrixD * ampsmods = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSM));
994 Int_t nrow = ampsmods->GetNrows();
995 Int_t ncol = ampsmods->GetNcols();
996 //printf("%s",Form(" ######## SM %i : nrow %i : ncol %i ##### \n", iSM, nrow, ncol));
997 for(Int_t row=0; row<nrow; row++) {
998 for(Int_t col=0; col<ncol; col++) {
999 amp = (*ampsmods)(row,col);
1003 if(keyPrint) AliDebug(2,Form("%s",Form(" ** nPhiSm %i : nEtaSM %i : row %2.2i : col %2.2i -> ", nPhiSM, nEtaSM, row, col)));
1004 if(nEtaSM == 0) { // positive Z
1005 jrow = 3*nPhiSM + row/nPhiModSum;
1006 jcol = 6 + col / nEtaModSum;
1007 } else { // negative Z
1008 if(iSM<=9) jrow = 3*nPhiSM + 2 - row/nPhiModSum;
1009 else jrow = 3*nPhiSM + 1 - row/nPhiModSum; // half size
1010 jcol = 5 - col / nEtaModSum;
1012 if(keyPrint) AliDebug(2,Form("%s",Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat) \n", jrow, jcol, amp)));
1014 (*jetMat)(jrow,jcol) += amp;
1015 ampSum += amp; // For controling
1016 } else if(amp<0.0) {
1017 AliDebug(1,Form(" jrow %2.2i : jcol %2.2i : amp %f (jetMat: amp<0) \n", jrow, jcol, amp));
1023 if(ampSum <= 0.0) AliDebug(1,Form("FillJetMatrixFromSMs","ampSum %f (<=0.0) ", ampSum));
1026 //____________________________________________________________________________
1027 void AliEMCALTrigger::MakeSlidingPatch(const TMatrixD &jm, const Int_t nPatchSize, TMatrixD &JetMax)
1029 // Sliding patch : nPatchSize x nPatchSize (OVERLAP)
1030 static int keyPrint = 0;
1031 if(keyPrint) AliDebug(2,Form(" AliEMCALTrigger::MakeSlidingPatch() was started \n"));
1032 Double_t ampCur = 0.0, e=0.0;
1033 ampJetMax(0,0) = 0.0;
1034 ampJetMax(3,0) = 0.0; // unused now
1035 ampJetMax(4,0) = ampJetMax(5,0) = 0.0;
1036 for(Int_t row=0; row<fAmpJetMatrix->GetNrows(); row ++) {
1037 for(Int_t col=0; col<fAmpJetMatrix->GetNcols(); col++) {
1039 // check on patch size
1040 if( (row+nPatchSize-1) < fAmpJetMatrix->GetNrows() && (col+nPatchSize-1) < fAmpJetMatrix->GetNcols()){
1041 for(Int_t i = 0 ; i < nPatchSize ; i++) {
1042 for(Int_t j = 0 ; j < nPatchSize ; j++) {
1043 ampCur += jm(row+i, col+j);
1045 } // end cycle on patch
1046 if(ampCur > ampJetMax(0,0)){
1047 ampJetMax(0,0) = ampCur;
1048 ampJetMax(1,0) = row;
1049 ampJetMax(2,0) = col;
1051 } // check on patch size
1054 if(keyPrint) AliDebug(2,Form(" ampJetMax %i row %2i->%2i col %2i->%2i \n", Int_t(ampJetMax(0,0)), Int_t(ampJetMax(1,0)), Int_t(ampJetMax(1,0))+nPatchSize-1, Int_t(ampJetMax(2,0)), Int_t(ampJetMax(2,0))+nPatchSize-1));
1056 Double_t eCorrJetMatrix=0.0;
1057 if(fVZER0Mult > 0.0) {
1058 // Correct patch energy (adc) and jet patch matrix energy
1059 Double_t meanAmpBG = GetMeanEmcalPatchEnergy(Int_t(fVZER0Mult), nPatchSize)/0.0153;
1060 ampJetMax(4,0) = ampJetMax(0,0);
1061 ampJetMax(5,0) = meanAmpBG;
1063 Double_t eCorr = ampJetMax(0,0) - meanAmpBG;
1064 AliDebug(2,Form(" ampJetMax(0,0) %f meanAmpBG %f eCorr %f : ampJetMax(4,0) %f \n",
1065 ampJetMax(0,0), meanAmpBG, eCorr, ampJetMax(5,0)));
1066 ampJetMax(0,0) = eCorr;
1068 eCorrJetMatrix = GetMeanEmcalEnergy(Int_t(fVZER0Mult)) / 208.;
1070 // Fill patch energy matrix
1071 for(int row=Int_t(ampJetMax(1,0)); row<Int_t(ampJetMax(1,0))+nPatchSize; row++) {
1072 for(int col=Int_t(ampJetMax(2,0)); col<Int_t(ampJetMax(2,0))+nPatchSize; col++) {
1073 e = Double_t(jm(row,col)*0.0153); // 0.0153 - hard coded now
1074 if(eCorrJetMatrix > 0.0) { // BG subtraction case
1075 e -= eCorrJetMatrix;
1076 fJetMatrixE->SetBinContent(row+1, col+1, e);
1077 } else if(e > 0.0) {
1078 fJetMatrixE->SetBinContent(row+1, col+1, e);
1082 // PrintJetMatrix();
1083 // Set the jet trigger(s), multiple threshold now, Nov 19,2007
1084 for(Int_t i=0; i<fNJetThreshold; i++ ) {
1085 if(ampJetMax(0,0) >= fL1JetThreshold[i]) {
1086 SetInput(GetNameOfJetTrigger(i));
1091 //____________________________________________________________________________
1092 Double_t AliEMCALTrigger::GetEmcalSumAmp() const
1094 // Return sum of amplidutes from EMCal
1095 // Used calibration coefficeint for transition to energy
1096 return fAmpJetMatrix >0 ?fAmpJetMatrix->Sum() :0.0;
1099 //____________________________________________________________________________
1100 void AliEMCALTrigger::PrintJetMatrix() const
1102 // fAmpJetMatrix : (17,12); // 17-phi(row), 12-eta(col)
1103 if(fAmpJetMatrix == 0) return;
1105 AliInfo(Form("\n #### jetMatrix : (%i,%i) ##### \n ",
1106 fAmpJetMatrix->GetNrows(), fAmpJetMatrix->GetNcols()));
1107 PrintMatrix(*fAmpJetMatrix);
1110 //____________________________________________________________________________
1111 void AliEMCALTrigger::PrintAmpTruMatrix(Int_t ind) const
1113 TMatrixD * tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1114 if(tru == 0) return;
1115 AliInfo(Form("\n #### Amp TRU matrix(%i) : (%i,%i) ##### \n ",
1116 ind, tru->GetNrows(), tru->GetNcols()));
1120 //____________________________________________________________________________
1121 void AliEMCALTrigger::PrintAmpSmMatrix(Int_t ind) const
1123 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(ind));
1125 AliInfo(Form("\n #### Amp SM matrix(%i) : (%i,%i) ##### \n ",
1126 ind, sm->GetNrows(), sm->GetNcols()));
1130 //____________________________________________________________________________
1131 void AliEMCALTrigger::PrintMatrix(const TMatrixD &mat) const
1133 for(Int_t col=0; col<mat.GetNcols(); col++) AliInfo(Form(" %3i ", col));
1134 AliInfo(Form("\n -- \n"));
1135 for(Int_t row=0; row<mat.GetNrows(); row++) {
1136 AliInfo(Form(" row:%2i ", row));
1137 for(Int_t col=0; col<mat.GetNcols(); col++) {
1138 AliInfo(Form(" %4i", (Int_t)mat(row,col)));
1144 //____________________________________________________________________________
1145 Bool_t AliEMCALTrigger::CheckConsistentOfMatrixes(const Int_t pri)
1147 Double_t sumSM = 0.0, smCur=0.0;
1148 Double_t sumTru=0.0, sumTruInSM = 0.0, truSum=0.0;
1149 // Bool_t key = kTRUE;
1151 for(Int_t i=0; i<fAmpSMods->GetEntries(); i++) {
1152 TMatrixD * sm = dynamic_cast<TMatrixD *>(fAmpSMods->At(i));
1158 for(Int_t itru=0; itru<3; itru++) { // Cycle on tru inside SM
1159 Int_t ind = 3*i + itru;
1160 TMatrixD *tru = dynamic_cast<TMatrixD *>(fAmpTrus->At(ind));
1162 truSum = tru->Sum();
1163 sumTruInSM += truSum;
1166 sumTru += sumTruInSM;
1168 if(sumTruInSM != smCur) {
1169 AliDebug(1,Form(" sm %i : smCur %f -> sumTruInSM %f \n", i, smCur, sumTruInSM));
1174 Double_t sumJetMat = fAmpJetMatrix->Sum();
1175 if(pri || sumSM != sumTru || sumSM != sumJetMat)
1176 AliDebug(1,Form(" sumSM %f : sumTru %f : sumJetMat %f \n", sumSM, sumTru, sumJetMat));
1177 if(sumSM != sumTru || sumSM != sumJetMat) return kFALSE;
1181 //____________________________________________________________________________
1182 void AliEMCALTrigger::Browse(TBrowser* b)
1184 if(&fInputs) b->Add(&fInputs);
1185 if(fAmpTrus) b->Add(fAmpTrus);
1186 if(fTimeRtrus) b->Add(fTimeRtrus);
1187 if(fAmpSMods) b->Add(fAmpSMods);
1188 if(fAmpJetMatrix) b->Add(fAmpJetMatrix);
1189 if(fJetMatrixE) b->Add(fJetMatrixE);