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