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