EMCAL/AliEMCALTrigger.cxx

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