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