EMCAL/AliEMCALTrigger.cxx

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