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