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