correction in C commentaries
[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 /* $Log$ */
17 /* Revision 1.7  2007/03/04 14:23:50  gustavo */
18 /* Trigger maximum amplitude patch value and postion stored in ESDs, possibility to isolate patch in SuperModule or TRU acceptance implemented */
19 /* */
20
21 //_________________________________________________________________________  
22 //
23 //  Class for trigger analysis.
24 //  Digits are grouped in TRU's  (Trigger Units). A TRU consists of 384 
25 //  cells ordered fNTRUPhi x fNTRUEta. The algorithm searches all possible 2x2 
26 //  and nxn (n is a multiple of 2) cell combinations per each TRU,  adding the 
27 //  digits amplitude and finding the maximum. If found, look if it is isolated.
28 //  Maxima are transformed in ADC time samples. Each time bin is compared to the trigger 
29 //  threshold until it is larger and then, triggers are set. Thresholds need to be fixed.  
30 //  Thresholds need to be fixed. Last 2 modules are half size in Phi, I considered 
31 //  that the number of TRU is maintained for the last modules but decision not taken. 
32 //  If different, then this must be changed. 
33 //  Usage:
34 //
35 //  //Inside the event loop
36 //  AliEMCALTrigger *tr = new AliEMCALTrigger();//Init Trigger
37 //  tr->SetL0Threshold(100); //Arbitrary threshold values
38 //  tr->SetL1JetLowPtThreshold(1000);
39 //  tr->SetL1JetMediumPtThreshold(10000);
40 //  tr->SetL1JetHighPtThreshold(20000);
41 //  ...
42 //  tr->Trigger(); //Execute Trigger
43 //  tr->Print(""); //Print results
44 //
45 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, CERN) 
46 //////////////////////////////////////////////////////////////////////////////
47
48
49 // --- ROOT system ---
50
51 // --- ALIROOT system ---
52 #include "AliRun.h" 
53 #include "AliRunLoader.h" 
54 #include "AliTriggerInput.h"
55 #include "AliEMCAL.h" 
56 #include "AliEMCALLoader.h" 
57 #include "AliEMCALDigit.h"
58 #include "AliEMCALTrigger.h" 
59 #include "AliEMCALGeometry.h"
60
61 ClassImp(AliEMCALTrigger)
62
63 //______________________________________________________________________
64 AliEMCALTrigger::AliEMCALTrigger()
65   : AliTriggerDetector(),
66     f2x2MaxAmp(-1), f2x2CellPhi(-1),  f2x2CellEta(-1),
67     f2x2SM(0),
68     fnxnMaxAmp(-1), fnxnCellPhi(-1),  fnxnCellEta(-1),
69     fnxnSM(0),
70     fADCValuesHighnxn(0),fADCValuesLownxn(0),
71     fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
72     fDigitsList(0),
73     fL0Threshold(100),fL1JetLowPtThreshold(200),
74     fL1JetMediumPtThreshold(500), fL1JetHighPtThreshold(1000),
75     fNTRU(3), fNTRUEta(3), fNTRUPhi(1), 
76     fNCellsPhi(24), fNCellsEta(16), 
77     fPatchSize(1),  fIsolPatchSize(1), 
78     f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1), 
79     f2x2AmpOutOfPatchThres(100000),  fnxnAmpOutOfPatchThres(100000), 
80     fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),  
81     fSimulation(kTRUE), fIsolateInSuperModule(kTRUE)
82 {
83   //ctor 
84   fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins];
85   fADCValuesLownxn  = 0x0; //new Int_t[fTimeBins];
86   fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins];
87   fADCValuesLow2x2  = 0x0; //new Int_t[fTimeBins];
88
89   SetName("EMCAL");
90   CreateInputs();
91    
92    //Print("") ; 
93 }
94
95
96
97 //____________________________________________________________________________
98 AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig) 
99   : AliTriggerDetector(trig),
100     f2x2MaxAmp(trig.f2x2MaxAmp), 
101     f2x2CellPhi(trig.f2x2CellPhi),  
102     f2x2CellEta(trig.f2x2CellEta),
103     f2x2SM(trig.f2x2SM),
104     fnxnMaxAmp(trig.fnxnMaxAmp), 
105     fnxnCellPhi(trig.fnxnCellPhi),  
106     fnxnCellEta(trig.fnxnCellEta),
107     fnxnSM(trig.fnxnSM),
108     fADCValuesHighnxn(trig.fADCValuesHighnxn),
109     fADCValuesLownxn(trig.fADCValuesLownxn),
110     fADCValuesHigh2x2(trig.fADCValuesHigh2x2),
111     fADCValuesLow2x2(trig.fADCValuesLow2x2),
112     fDigitsList(trig.fDigitsList),
113     fL0Threshold(trig.fL0Threshold),
114     fL1JetLowPtThreshold(trig.fL1JetLowPtThreshold),
115     fL1JetMediumPtThreshold(trig.fL1JetMediumPtThreshold), 
116     fL1JetHighPtThreshold(trig.fL1JetHighPtThreshold),
117     fNTRU(trig.fNTRU), fNTRUEta(trig.fNTRUEta), fNTRUPhi(trig.fNTRUPhi), 
118     fNCellsPhi(trig.fNCellsPhi), fNCellsEta(trig.fNCellsEta), 
119     fPatchSize(trig.fPatchSize),
120     fIsolPatchSize(trig.fIsolPatchSize), 
121     f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch), 
122     fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch), 
123     f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres),  
124     fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres), 
125     fIs2x2Isol(trig.fIs2x2Isol),
126     fIsnxnIsol(trig.fIsnxnIsol),  
127     fSimulation(trig.fSimulation),
128     fIsolateInSuperModule(trig.fIsolateInSuperModule)
129 {
130   // cpy ctor
131 }
132
133 //----------------------------------------------------------------------
134 void AliEMCALTrigger::CreateInputs()
135 {
136    // inputs 
137    
138    // Do not create inputs again!!
139    if( fInputs.GetEntriesFast() > 0 ) return;
140    
141    fInputs.AddLast( new AliTriggerInput( "EMCAL_L0",       "EMCAL L0", 0x02 ) );
142    fInputs.AddLast( new AliTriggerInput( "EMCAL_JetHPt_L1","EMCAL Jet High Pt L1",   0x04 ) );
143    fInputs.AddLast( new AliTriggerInput( "EMCAL_JetMPt_L1","EMCAL Jet Medium Pt L1", 0x08 ) );
144    fInputs.AddLast( new AliTriggerInput( "EMCAL_JetLPt_L1","EMCAL Jet Low Pt L1",    0x016 ) );
145  
146 }
147
148 //____________________________________________________________________________
149 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) {
150
151   //Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch, 
152   //inside isolation patch . iPatchType = 0 means calculation for 2x2 patch, 
153   //iPatchType = 1 means calculation for nxn patch.
154   //In the next table there is an example of the different options of patch size and isolation patch size:
155   //                                                                                 Patch Size (fPatchSize)
156   //                                                             0                          1                                  2
157   //          fIsolPatchSize                 2x2 (not overlap)   4x4 (overlapped)        6x6(overlapped) ...
158   //                   1                                       4x4                      8x8                              10x10
159   //                   2                                       6x6                     12x12                           14x14    
160   //                   3                                       8x8                     16x16                           18x18
161                           
162   Bool_t b = kFALSE;
163   Float_t amp = 0;
164  
165   //Get matrix of TRU or Module with maximum amplitude patch.
166   Int_t itru = mtru+iSM*fNTRU ; //number of tru, min 0 max 8*5.
167   TMatrixD * ampmatrix   = 0x0;
168   Int_t colborder = 0;
169   Int_t rowborder = 0;
170   
171   if(fIsolateInSuperModule){
172     ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
173     rowborder = fNCellsPhi*fNTRUPhi;
174     colborder = fNCellsEta*fNTRUEta;
175     AliDebug(2,"Isolate trigger in Module");
176   }
177   else{
178     ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
179     rowborder = fNCellsPhi;
180     colborder = fNCellsEta;
181     AliDebug(2,"Isolate trigger in TRU");
182   }
183   
184   //Define patch cells
185   Int_t isolcells   = fIsolPatchSize*(1+iPatchType);
186   Int_t ipatchcells = 2*(1+fPatchSize*iPatchType);
187   Int_t minrow      = maxphi - isolcells;
188   Int_t mincol      = maxeta - isolcells;
189   Int_t maxrow      = maxphi + isolcells + ipatchcells;
190   Int_t maxcol      = maxeta + isolcells + ipatchcells;
191
192   AliDebug(2,Form("Number of added Isol Cells %d, Patch Size %d",isolcells, ipatchcells));
193   AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
194   
195   if(minrow < 0 || mincol < 0 || maxrow > rowborder || maxcol > colborder){
196     AliDebug(1,Form("Out of Module/TRU range, cannot isolate patch"));
197     return kFALSE;
198   }
199   
200   //Add amplitudes in all isolation patch
201   for(Int_t irow = minrow ; irow <  maxrow; irow ++)
202     for(Int_t icol = mincol ; icol < maxcol ; icol ++)
203       amp += (*ampmatrix)(irow,icol);
204   
205   AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
206
207   if(amp < maxamp){
208     AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
209     return kFALSE;
210   }
211   else
212     amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
213   
214   AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
215
216   //Fill isolation amplitude data member and say if patch is isolated.
217   if(iPatchType == 0){ //2x2 case
218     f2x2AmpOutOfPatch = amp;   
219     if(amp < f2x2AmpOutOfPatchThres)
220       b=kTRUE;
221   }
222   else  if(iPatchType == 1){ //nxn case
223     fnxnAmpOutOfPatch = amp;   
224     if(amp < fnxnAmpOutOfPatchThres)
225       b=kTRUE;
226   }
227
228   return b;
229
230 }
231
232 //____________________________________________________________________________
233 void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD *ampmax2, TMatrixD *ampmaxn){
234   
235   //Sums energy of all possible 2x2 (L0) and nxn (L1) cells per each TRU. 
236   //Fast signal in the experiment is given by 2x2 cells, 
237   //for this reason we loop inside the TRU cells by 2. 
238
239   //Declare and initialize variables
240   Int_t nCellsPhi  = fNCellsPhi;//geom->GetNPhi()*2/geom->GetNTRUPhi() ;
241   if(isupermod > 9)
242     nCellsPhi =  nCellsPhi / 2 ; //Half size SM. Not Final.
243   // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
244   Int_t nCellsEta  = fNCellsEta ;//geom->GetNEta()*2/geom->GetNTRUEta() ;
245   // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
246   //Int_t nTRU          = geom->GeNTRU();//3 TRU per super module
247
248   Float_t amp2 = 0 ;
249   Float_t ampn = 0 ; 
250   for(Int_t i = 0; i < 4; i++){
251     for(Int_t j = 0; j < fNTRU; j++){
252       (*ampmax2)(i,j) = -1;
253       (*ampmaxn)(i,j) = -1;
254     }
255   }
256
257   //Create matrix that will contain 2x2 amplitude sums
258   //used to calculate the nxn sums
259   TMatrixD  * tru2x2 = new TMatrixD(nCellsPhi/2,nCellsEta/2) ;
260   for(Int_t i = 0; i < nCellsPhi/2; i++)
261     for(Int_t j = 0; j < nCellsEta/2; j++)
262       (*tru2x2)(i,j) = -1;
263   
264   //Loop over all TRUS in a supermodule
265   for(Int_t itru = 0 + isupermod * fNTRU ; itru < (isupermod+1)*fNTRU ; itru++) {
266     TMatrixD * amptru   = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
267     TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
268     Int_t mtru = itru-isupermod*fNTRU ; //Number of TRU in Supermodule
269    
270     //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
271     for(Int_t irow = 0 ; irow <  nCellsPhi; irow += 2){ 
272       for(Int_t icol = 0 ; icol < nCellsEta ; icol += 2){
273         amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
274           (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
275
276         //Fill matrix with added 2x2 cells for use in nxn sums
277         (*tru2x2)(irow/2,icol/2) = amp2 ;
278         //Select 2x2 maximum sums to select L0 
279         if(amp2 > (*ampmax2)(0,mtru)){
280           (*ampmax2)(0,mtru) = amp2 ; 
281           (*ampmax2)(1,mtru) = irow;
282           (*ampmax2)(2,mtru) = icol;
283         }
284       }
285     }
286     
287     //Find most recent time in the selected 2x2 cell
288     (*ampmax2)(3,mtru) = 1 ;
289     Int_t row2 =  static_cast <Int_t> ((*ampmax2)(1,mtru));
290     Int_t col2 =  static_cast <Int_t> ((*ampmax2)(2,mtru));
291     for(Int_t i = 0; i<2; i++){
292       for(Int_t j = 0; j<2; j++){
293         if((*amptru)(row2+i,col2+j) > 0 &&  (*timeRtru)(row2+i,col2+j)> 0){       
294           if((*timeRtru)(row2+i,col2+j) <  (*ampmax2)(3,mtru)  )
295             (*ampmax2)(3,mtru) =  (*timeRtru)(row2+i,col2+j);
296         }
297       }
298     }
299
300     //Sliding nxn, add nxn amplitudes  (OVERLAP)
301     if(fPatchSize > 0){
302       for(Int_t irow = 0 ; irow <  nCellsPhi/2; irow++){ 
303         for(Int_t icol = 0 ; icol < nCellsEta/2 ; icol++){
304           ampn = 0;
305           if( (irow+fPatchSize) < nCellsPhi/2 && (icol+fPatchSize) < nCellsEta/2){//Avoid exit the TRU
306             for(Int_t i = 0 ; i <= fPatchSize ; i++)
307               for(Int_t j = 0 ; j <= fPatchSize ; j++)
308                 ampn += (*tru2x2)(irow+i,icol+j);
309             //Select nxn maximum sums to select L1 
310             if(ampn > (*ampmaxn)(0,mtru)){
311               (*ampmaxn)(0,mtru) = ampn ; 
312               (*ampmaxn)(1,mtru) = irow*2;
313               (*ampmaxn)(2,mtru) = icol*2;
314             }
315           }
316         }
317       }
318       
319       //Find most recent time in selected nxn cell
320       (*ampmaxn)(3,mtru) = 1 ;
321       Int_t rown =  static_cast <Int_t> ((*ampmaxn)(1,mtru));
322       Int_t coln =  static_cast <Int_t> ((*ampmaxn)(2,mtru));
323       for(Int_t i = 0; i<4*fPatchSize; i++){
324         for(Int_t j = 0; j<4*fPatchSize; j++){
325           if( (rown+i) < nCellsPhi && (coln+j) < nCellsEta/2){//Avoid exit the TRU
326             if((*amptru)(rown+i,coln+j) > 0 &&  (*timeRtru)(rown+i,coln+j)> 0){
327               if((*timeRtru)(rown+i,coln+j) <  (*ampmaxn)(3,mtru)  )
328                 (*ampmaxn)(3,mtru) =  (*timeRtru)(rown+i,coln+j);
329             }
330           }
331         }
332       }
333     }
334     else {  
335         (*ampmaxn)(0,mtru) =  (*ampmax2)(0,mtru); 
336         (*ampmaxn)(1,mtru) =  (*ampmax2)(1,mtru);
337         (*ampmaxn)(2,mtru) =  (*ampmax2)(2,mtru);
338         (*ampmaxn)(3,mtru) =  (*ampmax2)(3,mtru);
339       }
340   }
341 }
342
343 //____________________________________________________________________________
344 void AliEMCALTrigger::Print(const Option_t * opt) const 
345 {
346   
347   //Prints main parameters
348   
349   if(! opt)
350     return;
351   AliTriggerInput* in = 0x0 ;
352
353   printf( "             Maximum Amplitude after Sliding Cell, \n") ; 
354   printf( "               -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
355           f2x2MaxAmp,f2x2SM) ; 
356   printf( "               -2x2 from row %d to row %d and from column %d to column %d\n", f2x2CellPhi, f2x2CellPhi+2, f2x2CellEta, f2x2CellEta+2) ; 
357   printf( "               -2x2 Isolation Patch %d x %d, Amplitude out of 2x2 patch is %f, threshold %f, Isolated? %d \n", 
358           2*fIsolPatchSize+2, 2*fIsolPatchSize+2,  f2x2AmpOutOfPatch,  f2x2AmpOutOfPatchThres,static_cast<Int_t> (fIs2x2Isol)) ; 
359   if(fPatchSize > 0){
360     printf( "             Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1));
361     printf( "               -nxn cells sum (overlapped)    : %10.2f, in Super Module %d\n",
362             fnxnMaxAmp,fnxnSM) ; 
363     printf( "               -nxn from row %d to row %d and from column %d to column %d\n", fnxnCellPhi, fnxnCellPhi+4*fPatchSize, fnxnCellEta, fnxnCellEta+4*fPatchSize) ; 
364     printf( "               -nxn Isolation Patch %d x %d, Amplitude out of nxn patch is %f, threshold %f, Isolated? %d \n", 
365             4*fIsolPatchSize+2*(fPatchSize+1),4*fIsolPatchSize+2*(fPatchSize+1) ,  fnxnAmpOutOfPatch,  fnxnAmpOutOfPatchThres,static_cast<Int_t> (fIsnxnIsol) ) ; 
366   }
367   
368   printf( "             Isolate in SuperModule? %d\n",  
369           fIsolateInSuperModule) ;  
370
371   printf( "             Threshold for LO %10.2f\n", 
372           fL0Threshold) ;  
373   in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
374   if(in->GetValue())
375     printf( "             *** EMCAL LO is set ***\n") ; 
376   
377   printf( "             Jet Low Pt Threshold for L1 %10.2f\n", 
378           fL1JetLowPtThreshold) ;
379   in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_JetLPt_L1" );
380   if(in->GetValue())
381     printf( "             *** EMCAL Jet Low Pt for L1 is set ***\n") ;
382
383   printf( "             Jet Medium Pt Threshold for L1 %10.2f\n", 
384           fL1JetMediumPtThreshold) ;  
385   in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_JetMPt_L1" );
386   if(in->GetValue())
387     printf( "             *** EMCAL Jet Medium Pt for L1 is set ***\n") ;
388
389   printf( "             Jet High Pt Threshold for L1 %10.2f\n", 
390           fL1JetHighPtThreshold) ;  
391   in = (AliTriggerInput*) fInputs.FindObject( "EMCAL_JetHPt_L1" );
392   if(in->GetValue())
393     printf( "             *** EMCAL Jet High Pt for L1 is set ***\n") ;
394
395 }
396
397 //____________________________________________________________________________
398 void AliEMCALTrigger::SetTriggers(const TClonesArray * ampmatrix,const Int_t iSM, 
399                                   const TMatrixD *ampmax2, 
400                                   const TMatrixD *ampmaxn, const AliEMCALGeometry *geom)  
401 {
402
403   //Checks the 2x2 and nxn maximum amplitude per each TRU and 
404   //compares with the different L0 and L1 triggers thresholds
405   Float_t max2[] = {-1,-1,-1,-1} ;
406   Float_t maxn[] = {-1,-1,-1,-1} ;
407   Int_t   mtru2  = -1 ;
408   Int_t   mtrun  = -1 ;
409
410   //Find maximum summed amplitude of all the TRU 
411   //in a Super Module
412     for(Int_t i = 0 ; i < fNTRU ; i++){
413       if(max2[0] < (*ampmax2)(0,i) ){
414         max2[0] =  (*ampmax2)(0,i) ; // 2x2 summed max amplitude
415         max2[1] =  (*ampmax2)(1,i) ; // corresponding phi position in TRU
416         max2[2] =  (*ampmax2)(2,i) ; // corresponding eta position in TRU
417         max2[3] =  (*ampmax2)(3,i) ; // corresponding most recent time
418         mtru2   = i ;
419       }
420       if(maxn[0] < (*ampmaxn)(0,i) ){
421         maxn[0] =  (*ampmaxn)(0,i) ; // nxn summed max amplitude
422         maxn[1] =  (*ampmaxn)(1,i) ; // corresponding phi position in TRU
423         maxn[2] =  (*ampmaxn)(2,i) ; // corresponding eta position in TRU
424         maxn[3] =  (*ampmaxn)(3,i) ; // corresponding most recent time
425         mtrun   = i ;
426       }
427     }
428
429   //--------Set max amplitude if larger than in other Super Modules------------
430   Float_t maxtimeR2 = -1 ;
431   Float_t maxtimeRn = -1 ;
432   AliRunLoader *rl  = AliRunLoader::GetRunLoader();
433   AliRun * gAlice   = rl->GetAliRun(); 
434   AliEMCAL * emcal  = (AliEMCAL*)gAlice->GetDetector("EMCAL");
435   Int_t nTimeBins = emcal->GetRawFormatTimeBins() ;
436
437   //Set max of 2x2 amplitudes and select L0 trigger
438   if(max2[0] > f2x2MaxAmp ){
439     f2x2MaxAmp  = max2[0] ;
440     f2x2SM      = iSM ;
441     maxtimeR2   = max2[3] ;
442     geom->GetCellPhiEtaIndexInSModuleFromTRUIndex(mtru2, 
443                                                   static_cast<Int_t>(max2[1]),
444                                                   static_cast<Int_t>(max2[2]),
445                                                   f2x2CellPhi,f2x2CellEta) ;
446     
447     //Isolated patch?
448     if(fIsolateInSuperModule)
449       fIs2x2Isol =  IsPatchIsolated(0, ampmatrix, iSM, mtru2,  f2x2MaxAmp, f2x2CellPhi,f2x2CellEta) ;
450     else
451       fIs2x2Isol =  IsPatchIsolated(0, ampmatrix, iSM, mtru2,  f2x2MaxAmp,  static_cast<Int_t>(max2[1]), static_cast<Int_t>(max2[2])) ;
452
453     //Transform digit amplitude in Raw Samples
454     fADCValuesLow2x2  = new Int_t[nTimeBins];
455     fADCValuesHigh2x2 = new Int_t[nTimeBins];
456     emcal->RawSampledResponse(maxtimeR2, f2x2MaxAmp, fADCValuesHigh2x2, fADCValuesLow2x2) ; 
457     
458     //Set Trigger Inputs, compare ADC time bins until threshold is attained
459     //Set L0
460     for(Int_t i = 0 ; i < nTimeBins ; i++){
461       if(fADCValuesHigh2x2[i] >= fL0Threshold          || fADCValuesLow2x2[i] >= fL0Threshold){
462         SetInput("EMCAL_L0") ;
463         break;
464       }
465     }
466   }
467   
468   //------------Set max of nxn amplitudes and select L1 trigger---------
469   if(maxn[0] > fnxnMaxAmp ){
470     fnxnMaxAmp  = maxn[0] ;
471     fnxnSM      = iSM ;
472     maxtimeRn   = maxn[3] ;
473     geom->GetCellPhiEtaIndexInSModuleFromTRUIndex(mtrun, 
474                                                   static_cast<Int_t>(maxn[1]),
475                                                   static_cast<Int_t>(maxn[2]),
476                                                   fnxnCellPhi,fnxnCellEta) ; 
477     
478     //Isolated patch?
479     if(fIsolateInSuperModule)
480       fIsnxnIsol =  IsPatchIsolated(1, ampmatrix, iSM, mtrun,  fnxnMaxAmp, fnxnCellPhi, fnxnCellEta) ;
481     else
482       fIsnxnIsol =  IsPatchIsolated(1, ampmatrix, iSM, mtrun,  fnxnMaxAmp,  static_cast<Int_t>(maxn[1]), static_cast<Int_t>(maxn[2])) ;
483     
484     //Transform digit amplitude in Raw Samples
485     fADCValuesHighnxn = new Int_t[nTimeBins];
486     fADCValuesLownxn  = new Int_t[nTimeBins];
487     emcal->RawSampledResponse(maxtimeRn, fnxnMaxAmp, fADCValuesHighnxn, fADCValuesLownxn) ;
488     
489     //Set Trigger Inputs, compare ADC time bins until threshold is attained
490     //SetL1 Low
491     for(Int_t i = 0 ; i < nTimeBins ; i++){
492       if(fADCValuesHighnxn[i] >= fL1JetLowPtThreshold  || fADCValuesLownxn[i] >= fL1JetLowPtThreshold){
493         SetInput("EMCAL_JetLPt_L1") ;
494         break; 
495       }
496     }
497     
498     //SetL1 Medium
499     for(Int_t i = 0 ; i < nTimeBins ; i++){
500       if(fADCValuesHighnxn[i] >= fL1JetMediumPtThreshold || fADCValuesLownxn[i] >= fL1JetMediumPtThreshold){
501         SetInput("EMCAL_JetMPt_L1") ;
502         break;
503       }
504     }
505     
506     //SetL1 High
507     for(Int_t i = 0 ; i < nTimeBins ; i++){
508       if(fADCValuesHighnxn[i] >= fL1JetHighPtThreshold || fADCValuesLownxn[i] >= fL1JetHighPtThreshold){
509         SetInput("EMCAL_JetHPt_L1") ;
510         break;
511       }
512     }
513   } 
514 }
515
516 //____________________________________________________________________________
517 void AliEMCALTrigger::Trigger() 
518 {
519   //Main Method to select triggers.
520   AliRunLoader *rl = AliRunLoader::GetRunLoader();
521   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
522     (rl->GetDetectorLoader("EMCAL"));
523   rl->LoadgAlice();  //Neede by calls to AliRun in SetTriggers
524  
525   //Load EMCAL Geometry
526   AliEMCALGeometry * geom =AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
527
528   if (geom==0)
529     AliFatal("Did not get geometry from EMCALLoader");
530   
531   //Define parameters
532   Int_t nSuperModules = geom->GetNumberOfSuperModules() ; //12 SM in EMCAL
533   fNTRU       = geom->GetNTRU();    //3 TRU per super module
534   fNTRUEta    = geom->GetNTRUEta(); //3 TRU in eta per super module
535   fNTRUPhi    = geom->GetNTRUPhi(); //1 TRU  in phi per super module
536   fNCellsPhi  = geom->GetNPhi()*2/geom->GetNTRUPhi() ;
537   fNCellsEta  = geom->GetNEta()*2/geom->GetNTRUEta() ;
538
539   //Intialize data members each time the trigger is called in event loop
540   f2x2MaxAmp = -1; f2x2CellPhi = -1;  f2x2CellEta = -1;
541   fnxnMaxAmp = -1; fnxnCellPhi = -1;  fnxnCellEta = -1;
542
543   //Take the digits list if simulation
544   if(fSimulation){
545     rl->LoadDigits("EMCAL");
546     fDigitsList = emcalLoader->Digits() ;
547   }
548   if(!fDigitsList)
549     AliFatal("Digits not found !") ;
550   
551   //Take the digits list 
552   
553   //Fill TRU Matrix  
554   TClonesArray * amptrus   = new TClonesArray("TMatrixD",1000);
555   TClonesArray * ampsmods  = new TClonesArray("TMatrixD",1000);
556   TClonesArray * timeRtrus = new TClonesArray("TMatrixD",1000);
557   
558   geom->FillTRU(fDigitsList, amptrus, ampsmods, timeRtrus) ;
559
560   //Do Cell Sliding and select Trigger
561   //Initialize varible that will contain maximum amplitudes and 
562   //its corresponding cell position in eta and phi, and time.
563   TMatrixD  * ampmax2 = new TMatrixD(4,fNTRU) ;
564   TMatrixD  * ampmaxn = new TMatrixD(4,fNTRU) ;
565   
566   for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
567     //Do 2x2 and nxn sums, select maximums. 
568     MakeSlidingCell(amptrus, timeRtrus, iSM, ampmax2, ampmaxn);
569     
570     //Set the trigger
571     if(fIsolateInSuperModule)
572       SetTriggers(ampsmods,iSM,ampmax2,ampmaxn,geom) ;
573     if(!fIsolateInSuperModule)
574       SetTriggers(amptrus,iSM,ampmax2,ampmaxn,geom) ;
575   }
576   
577   //Print();
578
579 }