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