EMCAL/AliEMCALTrigger.cxx

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