EMCAL/AliEMCALTrigger.cxx

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