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