Trigger maximum amplitude patch value and postion stored in ESDs, possibility to...
authorgustavo <gustavo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 4 Mar 2007 14:23:50 +0000 (14:23 +0000)
committergustavo <gustavo@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sun, 4 Mar 2007 14:23:50 +0000 (14:23 +0000)
EMCAL/AliEMCALGeometry.cxx
EMCAL/AliEMCALGeometry.h
EMCAL/AliEMCALReconstructor.cxx
EMCAL/AliEMCALTrigger.cxx
EMCAL/AliEMCALTrigger.h

index 7c130eb..dc1863c 100644 (file)
@@ -626,17 +626,15 @@ void AliEMCALGeometry::DefineSamplingFraction()
 }
 
 //____________________________________________________________________________
-void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * timeRmatrix) {
-
+void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampmatrix, TClonesArray * ampmatrixsmod, TClonesArray * timeRmatrix) {
 
 //  Orders digits ampitudes list in fNTRU TRUs (384 cells) per supermodule. 
 //  Each TRU is a TMatrixD, and they are kept in TClonesArrays. The number of 
 //  TRU in phi is fNTRUPhi, and the number of TRU in eta is fNTRUEta.
 //  Last 2 modules are half size in Phi, I considered that the number of TRU
 //  is maintained for the last modules but decision not taken. If different, 
-//  then this must be changed. 
+//  then this must be changed. Also fill a matrix with all amplitudes in supermodule for isolation studies. 
  
-
   //Check data members
 
   if(fNTRUEta*fNTRUPhi != fNTRU)
@@ -647,15 +645,16 @@ void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampma
   Int_t nCellsPhi  = fNPhi*2/fNTRUPhi;
   Int_t nCellsPhi2 = fNPhi/fNTRUPhi; //HalfSize modules
   Int_t nCellsEta  = fNZ*2/fNTRUEta;
-  Int_t id      = -1; 
-  Float_t amp   = -1;
-  Float_t timeR = -1;
-  Int_t iSupMod = -1;
+
+  Int_t id       = -1; 
+  Float_t amp    = -1;
+  Float_t timeR  = -1;
+  Int_t iSupMod  = -1;
   Int_t nModule  = -1;
-  Int_t nIphi   = -1;
-  Int_t nIeta   = -1;
-  Int_t iphi    = -1;
-  Int_t ieta    = -1;
+  Int_t nIphi    = -1;
+  Int_t nIeta    = -1;
+  Int_t iphi     = -1;
+  Int_t ieta     = -1;
 
   //List of TRU matrices initialized to 0.
   for(Int_t k = 0; k < fNTRU*fNumberOfSuperModules; k++){
@@ -671,6 +670,17 @@ void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampma
     new((*timeRmatrix)[k]) TMatrixD(*timeRtrus) ; 
   }
   
+  //List of Modules matrices initialized to 0.
+  for(Int_t k = 0; k < fNumberOfSuperModules ; k++){
+    TMatrixD  * ampsmods   = new TMatrixD( fNPhi*2, fNZ*2) ;
+    for(Int_t i = 0; i <  fNPhi*2; i++){
+      for(Int_t j = 0; j <  fNZ*2; j++){
+       (*ampsmods)(i,j)   = 0.0;
+      }
+    }
+    new((*ampmatrixsmod)[k])   TMatrixD(*ampsmods) ;
+  }
+
   AliEMCALDigit * dig ;
   
   //Digits loop to fill TRU matrices with amplitudes.
@@ -714,7 +724,11 @@ void AliEMCALGeometry::FillTRU(const TClonesArray * digits, TClonesArray * ampma
     
     (*amptrus)(irow,icol) = amp ;
     (*timeRtrus)(irow,icol) = timeR ;
-
+    
+    //####################SUPERMODULE MATRIX ##################
+    TMatrixD * ampsmods   = dynamic_cast<TMatrixD *>(ampmatrixsmod->At(iSupMod)) ;
+    (*ampsmods)(iphi,ieta)   = amp ;
+    
   }
 }
 
index 7742645..946dcbd 100644 (file)
@@ -49,7 +49,7 @@ public:
   virtual void Browse(TBrowser* b);
   virtual Bool_t  IsFolder() const;
 
-  void FillTRU(const TClonesArray * digits, TClonesArray * amptru, TClonesArray * timeRtru)  ; //Fills Trigger Unit matrices with digit amplitudes and time
+  void FillTRU(const TClonesArray * digits, TClonesArray * amptru, TClonesArray * ampmod, TClonesArray * timeRtru)  ; //Fills Trigger Unit matrices with digit amplitudes and time
   void GetCellPhiEtaIndexInSModuleFromTRUIndex(Int_t itru, Int_t iphitru, Int_t ietatru, Int_t &ietaSM, Int_t &iphiSM) const ; // Tranforms Eta-Phi Cell index in TRU into Eta-Phi index in Super Module
   
   // Have to call GetTransformationForSM() before calculation global charachteristics 
index 442fcbe..eb1e6f0 100644 (file)
@@ -34,6 +34,7 @@
 #include "AliEMCALClusterizerv1.h"
 #include "AliEMCALRecPoint.h"
 #include "AliEMCALPID.h"
+#include "AliEMCALTrigger.h"
 #include "AliRawReader.h"
 
 
@@ -125,6 +126,58 @@ void AliEMCALReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
   //  Int_t nRP=0, nPC=0; // in input
   esd->SetFirstEMCALCluster(esd->GetNumberOfCaloClusters()); // Put after Phos clusters 
 
+  //#########Calculate trigger and set trigger info###########
+  AliEMCALTrigger tr ;
+  //   tr.SetPatchSize(1);//create 4x4 patches
+  tr.Trigger();
+  
+  Float_t maxAmp2x2  = tr.Get2x2MaxAmplitude();
+  Float_t maxAmpnxn  = tr.GetnxnMaxAmplitude();
+  Float_t ampOutOfPatch2x2  = tr.Get2x2AmpOutOfPatch() ;
+  Float_t ampOutOfPatchnxn  = tr.GetnxnAmpOutOfPatch() ;
+
+  AliEMCALGeometry * geom =  AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
+
+  Int_t iSM2x2      = tr.Get2x2SuperModule();
+  Int_t iSMnxn      = tr.GetnxnSuperModule();
+  Int_t iCellPhi2x2 = tr.Get2x2CellPhi();
+  Int_t iCellPhinxn = tr.GetnxnCellPhi();
+  Int_t iCellEta2x2 = tr.Get2x2CellEta();
+  Int_t iCellEtanxn = tr.GetnxnCellEta();
+
+  AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",  maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCellPhi2x2, iCellEta2x2));
+  AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",  maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCellPhinxn, iCellEtanxn));
+
+  TVector3    pos2x2(-1,-1,-1);
+  TVector3    posnxn(-1,-1,-1);
+
+  Int_t iAbsId2x2 = geom->GetAbsCellIdFromCellIndexes( iSM2x2, iCellPhi2x2, iCellEta2x2) ;
+  Int_t iAbsIdnxn = geom->GetAbsCellIdFromCellIndexes( iSMnxn, iCellPhinxn, iCellEtanxn) ;
+  geom->GetGlobal(iAbsId2x2, pos2x2);
+  geom->GetGlobal(iAbsIdnxn, posnxn);
+  
+  TArrayF triggerPosition(6);
+  triggerPosition[0] = pos2x2(0) ;   
+  triggerPosition[1] = pos2x2(1) ;   
+  triggerPosition[2] = pos2x2(2) ;  
+  triggerPosition[3] = posnxn(0) ;   
+  triggerPosition[4] = posnxn(1) ;   
+  triggerPosition[5] = posnxn(2) ;  
+
+  TArrayF triggerAmplitudes(4);
+  triggerAmplitudes[0] = maxAmp2x2 ;   
+  triggerAmplitudes[1] = ampOutOfPatch2x2 ;    
+  triggerAmplitudes[2] = maxAmpnxn ;   
+  triggerAmplitudes[3] = ampOutOfPatchnxn ;   
+
+  esd->AddEMCALTriggerPosition(triggerPosition);
+  esd->AddEMCALTriggerAmplitudes(triggerAmplitudes);
+  
+  //######################################
+
+  //Fill CaloClusters 
+
   for (Int_t iClust = 0 ; iClust < nClusters ; iClust++) {
     const AliEMCALRecPoint * clust = emcalLoader->RecPoint(iClust);
     //if(clust->GetClusterType()== AliESDCaloCluster::kClusterv1) nRP++; else nPC++;
index 257a8d0..339a1cb 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 /* $Id$ */
-/* $Log $ */
+/* $Log$ */
 
 //_________________________________________________________________________  
 //
 //  Class for trigger analysis.
-//  Digits are grouped in TRU's (384 cells ordered fNTRUPhi x fNTRUEta). 
-//  The algorithm searches all possible 2x2 and nxn (n is a multiple of 4) cell 
-//  combinations per each TRU,  adding the digits amplitude and finding the 
-//  maximum. Maxima are compared to triggers threshold and they are set. 
+//  Digits are grouped in TRU's  (Trigger Units). A TRU consists of 384 
+//  cells ordered fNTRUPhi x fNTRUEta. The algorithm searches all possible 2x2 
+//  and nxn (n is a multiple of 2) cell combinations per each TRU,  adding the 
+//  digits amplitude and finding the maximum. If found, look if it is isolated.
+//  Maxima are transformed in ADC time samples. Each time bin is compared to the trigger 
+//  threshold until it is larger and then, triggers are set. Thresholds need to be fixed.  
 //  Thresholds need to be fixed. Last 2 modules are half size in Phi, I considered 
 //  that the number of TRU is maintained for the last modules but decision not taken. 
 //  If different, then this must be changed. 
@@ -33,6 +35,7 @@
 //  tr->SetL1JetLowPtThreshold(1000);
 //  tr->SetL1JetMediumPtThreshold(10000);
 //  tr->SetL1JetHighPtThreshold(20000);
+//  ...
 //  tr->Trigger(); //Execute Trigger
 //  tr->Print(""); //Print results
 //
 
 
 // --- ROOT system ---
-//#include "TMatrixD.h"
 
 // --- ALIROOT system ---
-
 #include "AliRun.h" 
 #include "AliRunLoader.h" 
 #include "AliTriggerInput.h"
@@ -63,14 +64,24 @@ AliEMCALTrigger::AliEMCALTrigger()
     f2x2SM(0),
     fnxnMaxAmp(-1), fnxnCellPhi(-1),  fnxnCellEta(-1),
     fnxnSM(0),
-    fADCValuesHighnxn(0x0),fADCValuesLownxn(0x0),
-    fADCValuesHigh2x2(0x0),fADCValuesLow2x2(0x0),
-    fDigitsList(0x0),
+    fADCValuesHighnxn(0),fADCValuesLownxn(0),
+    fADCValuesHigh2x2(0),fADCValuesLow2x2(0),
+    fDigitsList(0),
     fL0Threshold(100),fL1JetLowPtThreshold(200),
     fL1JetMediumPtThreshold(500), fL1JetHighPtThreshold(1000),
-    fPatchSize(1), fSimulation(kTRUE)
+    fNTRU(3), fNTRUEta(3), fNTRUPhi(1), 
+    fNCellsPhi(24), fNCellsEta(16), 
+    fPatchSize(1),  fIsolPatchSize(1), 
+    f2x2AmpOutOfPatch(-1), fnxnAmpOutOfPatch(-1), 
+    f2x2AmpOutOfPatchThres(100000),  fnxnAmpOutOfPatchThres(100000), 
+    fIs2x2Isol(kFALSE), fIsnxnIsol(kFALSE),  
+    fSimulation(kTRUE), fIsolateInSuperModule(kTRUE)
 {
   //ctor 
+  fADCValuesHighnxn = 0x0; //new Int_t[fTimeBins];
+  fADCValuesLownxn  = 0x0; //new Int_t[fTimeBins];
+  fADCValuesHigh2x2 = 0x0; //new Int_t[fTimeBins];
+  fADCValuesLow2x2  = 0x0; //new Int_t[fTimeBins];
 
   SetName("EMCAL");
   CreateInputs();
@@ -100,7 +111,18 @@ AliEMCALTrigger::AliEMCALTrigger(const AliEMCALTrigger & trig)
     fL1JetLowPtThreshold(trig.fL1JetLowPtThreshold),
     fL1JetMediumPtThreshold(trig.fL1JetMediumPtThreshold), 
     fL1JetHighPtThreshold(trig.fL1JetHighPtThreshold),
-    fPatchSize(trig.fPatchSize), fSimulation(trig.fSimulation)
+    fNTRU(trig.fNTRU), fNTRUEta(trig.fNTRUEta), fNTRUPhi(trig.fNTRUPhi), 
+    fNCellsPhi(trig.fNCellsPhi), fNCellsEta(trig.fNCellsEta), 
+    fPatchSize(trig.fPatchSize),
+    fIsolPatchSize(trig.fIsolPatchSize), 
+    f2x2AmpOutOfPatch(trig.f2x2AmpOutOfPatch), 
+    fnxnAmpOutOfPatch(trig.fnxnAmpOutOfPatch), 
+    f2x2AmpOutOfPatchThres(trig.f2x2AmpOutOfPatchThres),  
+    fnxnAmpOutOfPatchThres(trig.fnxnAmpOutOfPatchThres), 
+    fIs2x2Isol(trig.fIs2x2Isol),
+    fIsnxnIsol(trig.fIsnxnIsol),  
+    fSimulation(trig.fSimulation),
+    fIsolateInSuperModule(trig.fIsolateInSuperModule)
 {
   // cpy ctor
 }
@@ -121,25 +143,109 @@ void AliEMCALTrigger::CreateInputs()
 }
 
 //____________________________________________________________________________
-void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD *ampmax2, TMatrixD *ampmaxn, AliEMCALGeometry *geom){
+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) {
+
+  //Calculate if the maximum patch found is isolated, find amplitude around maximum (2x2 or nxn) patch, 
+  //inside isolation patch . iPatchType = 0 means calculation for 2x2 patch, 
+  //iPatchType = 1 means calculation for nxn patch.
+  //In the next table there is an example of the different options of patch size and isolation patch size:
+  //                                                                                 Patch Size (fPatchSize)
+  //                                                             0                          1                                  2
+  //          fIsolPatchSize                 2x2 (not overlap)   4x4 (overlapped)        6x6(overlapped) ...
+  //                   1                                       4x4                      8x8                              10x10
+  //                   2                                       6x6                     12x12                           14x14    
+  //                   3                                       8x8                     16x16                           18x18
+                          
+  Bool_t b = kFALSE;
+  Float_t amp = 0;
+  //Get matrix of TRU or Module with maximum amplitude patch.
+  Int_t itru = mtru+iSM*fNTRU ; //number of tru, min 0 max 8*5.
+  TMatrixD * ampmatrix   = 0x0;
+  Int_t colborder = 0;
+  Int_t rowborder = 0;
+  
+  if(fIsolateInSuperModule){
+    ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(iSM)) ;
+    rowborder = fNCellsPhi*fNTRUPhi;
+    colborder = fNCellsEta*fNTRUEta;
+    AliDebug(2,"Isolate trigger in Module");
+  }
+  else{
+    ampmatrix = dynamic_cast<TMatrixD *>(ampmatrixes->At(itru)) ;
+    rowborder = fNCellsPhi;
+    colborder = fNCellsEta;
+    AliDebug(2,"Isolate trigger in TRU");
+  }
+  
+  //Define patch cells
+  Int_t isolcells   = fIsolPatchSize*(1+iPatchType);
+  Int_t ipatchcells = 2*(1+fPatchSize*iPatchType);
+  Int_t minrow      = maxphi - isolcells;
+  Int_t mincol      = maxeta - isolcells;
+  Int_t maxrow      = maxphi + isolcells + ipatchcells;
+  Int_t maxcol      = maxeta + isolcells + ipatchcells;
+
+  AliDebug(2,Form("Number of added Isol Cells %d, Patch Size %d",isolcells, ipatchcells));
+  AliDebug(2,Form("Patch: minrow %d, maxrow %d, mincol %d, maxcol %d",minrow,maxrow,mincol,maxcol));
+  
+  if(minrow < 0 || mincol < 0 || maxrow > rowborder || maxcol > colborder){
+    AliDebug(1,Form("Out of Module/TRU range, cannot isolate patch"));
+    return kFALSE;
+  }
+  
+  //Add amplitudes in all isolation patch
+  for(Int_t irow = minrow ; irow <  maxrow; irow ++)
+    for(Int_t icol = mincol ; icol < maxcol ; icol ++)
+      amp += (*ampmatrix)(irow,icol);
+  
+  AliDebug(2,Form("Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
+
+  if(amp < maxamp){
+    AliError(Form("Bad sum: Type %d, Maximum amplitude %f, patch+isol square %f",iPatchType, maxamp, amp));
+    return kFALSE;
+  }
+  else
+    amp-=maxamp; //Calculate energy in isolation patch that do not comes from maximum patch.
+  
+  AliDebug(2, Form("Maximum amplitude %f, Out of patch %f",maxamp, amp));
+
+  //Fill isolation amplitude data member and say if patch is isolated.
+  if(iPatchType == 0){ //2x2 case
+    f2x2AmpOutOfPatch = amp;   
+    if(amp < f2x2AmpOutOfPatchThres)
+      b=kTRUE;
+  }
+  else  if(iPatchType == 1){ //nxn case
+    fnxnAmpOutOfPatch = amp;   
+    if(amp < fnxnAmpOutOfPatchThres)
+      b=kTRUE;
+  }
+
+  return b;
+
+}
+
+//____________________________________________________________________________
+void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus, const Int_t isupermod,TMatrixD *ampmax2, TMatrixD *ampmaxn){
   
   //Sums energy of all possible 2x2 (L0) and nxn (L1) cells per each TRU. 
   //Fast signal in the experiment is given by 2x2 cells, 
   //for this reason we loop inside the TRU cells by 2. 
 
   //Declare and initialize variables
-  Int_t nCellsPhi  = geom->GetNPhi()*2/geom->GetNTRUPhi() ;
+  Int_t nCellsPhi  = fNCellsPhi;//geom->GetNPhi()*2/geom->GetNTRUPhi() ;
   if(isupermod > 9)
     nCellsPhi =  nCellsPhi / 2 ; //Half size SM. Not Final.
   // 12(tow)*2(cell)/1 TRU, cells in Phi in one TRU
-  Int_t nCellsEta  = geom->GetNEta()*2/geom->GetNTRUEta() ;
+  Int_t nCellsEta  = fNCellsEta ;//geom->GetNEta()*2/geom->GetNTRUEta() ;
   // 24(mod)*2(tower)/3 TRU, cells in Eta in one TRU
-  Int_t nTRU          = geom->GetNTRU();//3 TRU per super module
+  //Int_t nTRU          = geom->GeNTRU();//3 TRU per super module
 
   Float_t amp2 = 0 ;
   Float_t ampn = 0 ; 
   for(Int_t i = 0; i < 4; i++){
-    for(Int_t j = 0; j < nTRU; j++){
+    for(Int_t j = 0; j < fNTRU; j++){
       (*ampmax2)(i,j) = -1;
       (*ampmaxn)(i,j) = -1;
     }
@@ -153,10 +259,10 @@ void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClone
       (*tru2x2)(i,j) = -1;
   
   //Loop over all TRUS in a supermodule
-  for(Int_t itru = 0 + isupermod * nTRU ; itru < (isupermod+1)*nTRU ; itru++) {
+  for(Int_t itru = 0 + isupermod * fNTRU ; itru < (isupermod+1)*fNTRU ; itru++) {
     TMatrixD * amptru   = dynamic_cast<TMatrixD *>(amptrus->At(itru)) ;
     TMatrixD * timeRtru = dynamic_cast<TMatrixD *>(timeRtrus->At(itru)) ;
-    Int_t mtru = itru-isupermod*nTRU ; //Number of TRU in Supermodule
+    Int_t mtru = itru-isupermod*fNTRU ; //Number of TRU in Supermodule
    
     //Sliding 2x2, add 2x2 amplitudes (NOT OVERLAP)
     for(Int_t irow = 0 ; irow <  nCellsPhi; irow += 2){ 
@@ -164,7 +270,7 @@ void AliEMCALTrigger::MakeSlidingCell(const TClonesArray * amptrus, const TClone
        amp2 = (*amptru)(irow,icol)+(*amptru)(irow+1,icol)+
          (*amptru)(irow,icol+1)+(*amptru)(irow+1,icol+1);
 
-       //Fill matrix with added 2x2 crystals for use in nxn sums
+       //Fill matrix with added 2x2 cells for use in nxn sums
        (*tru2x2)(irow/2,icol/2) = amp2 ;
        //Select 2x2 maximum sums to select L0 
        if(amp2 > (*ampmax2)(0,mtru)){
@@ -245,14 +351,20 @@ void AliEMCALTrigger::Print(const Option_t * opt) const
   printf( "               -2x2 cells sum (not overlapped): %10.2f, in Super Module %d\n",
          f2x2MaxAmp,f2x2SM) ; 
   printf( "               -2x2 from row %d to row %d and from column %d to column %d\n", f2x2CellPhi, f2x2CellPhi+2, f2x2CellEta, f2x2CellEta+2) ; 
+  printf( "               -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)) ; 
   if(fPatchSize > 0){
-    printf( "             Patch Size, n x n: %d x %d cells\n",4*fPatchSize, 4*fPatchSize);
+    printf( "             Patch Size, n x n: %d x %d cells\n",2*(fPatchSize+1), 2*(fPatchSize+1));
     printf( "               -nxn cells sum (overlapped)    : %10.2f, in Super Module %d\n",
            fnxnMaxAmp,fnxnSM) ; 
     printf( "               -nxn from row %d to row %d and from column %d to column %d\n", fnxnCellPhi, fnxnCellPhi+4*fPatchSize, fnxnCellEta, fnxnCellEta+4*fPatchSize) ; 
+    printf( "               -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) ) ; 
   }
+  
+  printf( "             Isolate in SuperModule? %d\n",  
+          fIsolateInSuperModule) ;  
+
   printf( "             Threshold for LO %10.2f\n", 
          fL0Threshold) ;  
   in = (AliTriggerInput*)fInputs.FindObject( "EMCAL_L0" );
@@ -280,33 +392,34 @@ void AliEMCALTrigger::Print(const Option_t * opt) const
 }
 
 //____________________________________________________________________________
-void AliEMCALTrigger::SetTriggers(const Int_t iSM, const TMatrixD *ampmax2, 
-                                 const TMatrixD *ampmaxn, AliEMCALGeometry *geom)  
+void AliEMCALTrigger::SetTriggers(const TClonesArray * ampmatrix,const Int_t iSM, 
+                                 const TMatrixD *ampmax2, 
+                                 const TMatrixD *ampmaxn, const AliEMCALGeometry *geom)  
 {
 
   //Checks the 2x2 and nxn maximum amplitude per each TRU and 
   //compares with the different L0 and L1 triggers thresholds
   Float_t max2[] = {-1,-1,-1,-1} ;
   Float_t maxn[] = {-1,-1,-1,-1} ;
-  Int_t   itru2  = -1 ;
-  Int_t   itrun  = -1 ;
+  Int_t   mtru2  = -1 ;
+  Int_t   mtrun  = -1 ;
 
   //Find maximum summed amplitude of all the TRU 
   //in a Super Module
-    for(Int_t i = 0 ; i < geom->GetNTRU() ; i++){
+    for(Int_t i = 0 ; i < fNTRU ; i++){
       if(max2[0] < (*ampmax2)(0,i) ){
        max2[0] =  (*ampmax2)(0,i) ; // 2x2 summed max amplitude
        max2[1] =  (*ampmax2)(1,i) ; // corresponding phi position in TRU
        max2[2] =  (*ampmax2)(2,i) ; // corresponding eta position in TRU
        max2[3] =  (*ampmax2)(3,i) ; // corresponding most recent time
-       itru2   = i ;
+       mtru2   = i ;
       }
       if(maxn[0] < (*ampmaxn)(0,i) ){
        maxn[0] =  (*ampmaxn)(0,i) ; // nxn summed max amplitude
        maxn[1] =  (*ampmaxn)(1,i) ; // corresponding phi position in TRU
        maxn[2] =  (*ampmaxn)(2,i) ; // corresponding eta position in TRU
        maxn[3] =  (*ampmaxn)(3,i) ; // corresponding most recent time
-       itrun   = i ;
+       mtrun   = i ;
       }
     }
 
@@ -323,11 +436,17 @@ void AliEMCALTrigger::SetTriggers(const Int_t iSM, const TMatrixD *ampmax2,
     f2x2MaxAmp  = max2[0] ;
     f2x2SM      = iSM ;
     maxtimeR2   = max2[3] ;
-    geom->GetCellPhiEtaIndexInSModuleFromTRUIndex(itru2, 
+    geom->GetCellPhiEtaIndexInSModuleFromTRUIndex(mtru2, 
                                                  static_cast<Int_t>(max2[1]),
                                                  static_cast<Int_t>(max2[2]),
                                                  f2x2CellPhi,f2x2CellEta) ;
     
+    //Isolated patch?
+    if(fIsolateInSuperModule)
+      fIs2x2Isol =  IsPatchIsolated(0, ampmatrix, iSM, mtru2,  f2x2MaxAmp, f2x2CellPhi,f2x2CellEta) ;
+    else
+      fIs2x2Isol =  IsPatchIsolated(0, ampmatrix, iSM, mtru2,  f2x2MaxAmp,  static_cast<Int_t>(max2[1]), static_cast<Int_t>(max2[2])) ;
+
     //Transform digit amplitude in Raw Samples
     fADCValuesLow2x2  = new Int_t[nTimeBins];
     fADCValuesHigh2x2 = new Int_t[nTimeBins];
@@ -348,10 +467,17 @@ void AliEMCALTrigger::SetTriggers(const Int_t iSM, const TMatrixD *ampmax2,
     fnxnMaxAmp  = maxn[0] ;
     fnxnSM      = iSM ;
     maxtimeRn   = maxn[3] ;
-    geom->GetCellPhiEtaIndexInSModuleFromTRUIndex(itrun, 
+    geom->GetCellPhiEtaIndexInSModuleFromTRUIndex(mtrun, 
                                                  static_cast<Int_t>(maxn[1]),
                                                  static_cast<Int_t>(maxn[2]),
                                                  fnxnCellPhi,fnxnCellEta) ; 
+    
+    //Isolated patch?
+    if(fIsolateInSuperModule)
+      fIsnxnIsol =  IsPatchIsolated(1, ampmatrix, iSM, mtrun,  fnxnMaxAmp, fnxnCellPhi, fnxnCellEta) ;
+    else
+      fIsnxnIsol =  IsPatchIsolated(1, ampmatrix, iSM, mtrun,  fnxnMaxAmp,  static_cast<Int_t>(maxn[1]), static_cast<Int_t>(maxn[2])) ;
+    
     //Transform digit amplitude in Raw Samples
     fADCValuesHighnxn = new Int_t[nTimeBins];
     fADCValuesLownxn  = new Int_t[nTimeBins];
@@ -391,20 +517,21 @@ void AliEMCALTrigger::Trigger()
   AliRunLoader *rl = AliRunLoader::GetRunLoader();
   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
     (rl->GetDetectorLoader("EMCAL"));
-
+  rl->LoadgAlice();  //Neede by calls to AliRun in SetTriggers
   //Load EMCAL Geometry
-  rl->LoadgAlice(); 
-  AliRun * gAlice = rl->GetAliRun(); 
-  AliEMCAL * emcal  = (AliEMCAL*)gAlice->GetDetector("EMCAL");
-  AliEMCALGeometry * geom = emcal->GetGeometry();
-  
+  AliEMCALGeometry * geom =AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
+
   if (geom==0)
     AliFatal("Did not get geometry from EMCALLoader");
-
-
+  
   //Define parameters
   Int_t nSuperModules = geom->GetNumberOfSuperModules() ; //12 SM in EMCAL
-  Int_t nTRU          = geom->GetNTRU();//3 TRU per super module
+  fNTRU       = geom->GetNTRU();    //3 TRU per super module
+  fNTRUEta    = geom->GetNTRUEta(); //3 TRU in eta per super module
+  fNTRUPhi    = geom->GetNTRUPhi(); //1 TRU  in phi per super module
+  fNCellsPhi  = geom->GetNPhi()*2/geom->GetNTRUPhi() ;
+  fNCellsEta  = geom->GetNEta()*2/geom->GetNTRUEta() ;
 
   //Intialize data members each time the trigger is called in event loop
   f2x2MaxAmp = -1; f2x2CellPhi = -1;  f2x2CellEta = -1;
@@ -422,21 +549,28 @@ void AliEMCALTrigger::Trigger()
   
   //Fill TRU Matrix  
   TClonesArray * amptrus   = new TClonesArray("TMatrixD",1000);
+  TClonesArray * ampsmods  = new TClonesArray("TMatrixD",1000);
   TClonesArray * timeRtrus = new TClonesArray("TMatrixD",1000);
-
-  geom->FillTRU(fDigitsList, amptrus, timeRtrus) ;
+  
+  geom->FillTRU(fDigitsList, amptrus, ampsmods, timeRtrus) ;
 
   //Do Cell Sliding and select Trigger
   //Initialize varible that will contain maximum amplitudes and 
   //its corresponding cell position in eta and phi, and time.
-  TMatrixD  * ampmax2 = new TMatrixD(4,nTRU) ;
-  TMatrixD  * ampmaxn = new TMatrixD(4,nTRU) ;
-
+  TMatrixD  * ampmax2 = new TMatrixD(4,fNTRU) ;
+  TMatrixD  * ampmaxn = new TMatrixD(4,fNTRU) ;
+  
   for(Int_t iSM = 0 ; iSM < nSuperModules ; iSM++) {
     //Do 2x2 and nxn sums, select maximums. 
-    MakeSlidingCell(amptrus, timeRtrus, iSM, ampmax2, ampmaxn, geom);
-  
+    MakeSlidingCell(amptrus, timeRtrus, iSM, ampmax2, ampmaxn);
+    
     //Set the trigger
-    SetTriggers(iSM, ampmax2, ampmaxn, geom) ;
+    if(fIsolateInSuperModule)
+      SetTriggers(ampsmods,iSM,ampmax2,ampmaxn,geom) ;
+    if(!fIsolateInSuperModule)
+      SetTriggers(amptrus,iSM,ampmax2,ampmaxn,geom) ;
   }
+  
+  //Print();
+
 }
index 2474e6d..c7b1d71 100644 (file)
@@ -3,16 +3,20 @@
 /* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
 
-/* $Id $ */
-/* $Log $ */
+/* $Id$ */
+/* $Log$ */
+
 //___________________________________________________________
 //  Class for trigger analysis.
+//
+//  -- Author: Gustavo Conesa & Yves Schutz (IFIC, SUBATECH, CERN)
 //  Digits are grouped in TRU's (Trigger Units). A TRU consist of 384 cells 
-//  ordered fNTRUPhi x fNTRUZ. The algorithm searches all possible 2x2 and 
-//  nxn (n multiple of 4) crystal combinations per each TRU, adding the digits 
-//  amplitude and finding the maximum. Maximums are transformed in adc time 
-//  samples. Each time bin is compared to the trigger threshold until it is larger 
-//  and then, triggers are set. Thresholds need to be fixed. 
+//  ordered fNTRUPhi x fNTRUEta matrix. The algorithm searches all possible 
+//  2x2 and nxn (n multiple of 4) crystal combinations per each TRU, adding the 
+//  digits amplitude and finding the maximum.  It is found is maximum is isolated. 
+//  Maxima are transformed in adc time samples. Each time bin is compared to the 
+//  trigger threshold until it is larger and then, triggers are set. 
+//  Thresholds need to be fixed. 
 //  Last 2 modules are half size in Phi, I considered that the number 
 //  of TRU is maintained for the last modules but final decision has not 
 //  been taken. If different, then this must to be changed. 
@@ -24,8 +28,9 @@
 //  tr->SetL1JetLowPtThreshold(1000);
 //  tr->SetL1JetMediumPtThreshold(10000);
 //  tr->SetL1JetHighPtThreshold(20000);
+//  ....
 //  tr->Trigger(); //Execute Trigger
-//  tr->Print(""); //Print results
+//  tr->Print("");  //Print data members after calculation.
 //
 //*-- Author: Gustavo Conesa & Yves Schutz (IFIC, SUBATECH, CERN)
      
@@ -46,6 +51,8 @@ class AliEMCALTrigger : public AliTriggerDetector {
   AliEMCALTrigger() ; //  ctor
   AliEMCALTrigger(const AliEMCALTrigger & trig) ; // cpy ctor
   virtual ~AliEMCALTrigger() {}; //virtual dtor
+
+
   virtual void    CreateInputs(); //Define trigger inputs for Central Trigger Processor
   void            Print(const Option_t * opt ="") const ;  
   virtual void    Trigger();  //Make EMCAL trigger
@@ -54,28 +61,39 @@ class AliEMCALTrigger : public AliTriggerDetector {
   const AliEMCALTrigger & operator = (const AliEMCALTrigger & ) {return *this;}
 
   //Getters
-  Float_t  Get2x2MaxAmplitude()  const {return f2x2MaxAmp ; }
-  Float_t  GetnxnMaxAmplitude()  const {return fnxnMaxAmp ; }
-  Int_t    Get2x2CellPhi()       const {return f2x2CellPhi ; }
-  Int_t    GetnxnCellPhi()       const {return fnxnCellPhi ; }
-  Int_t    Get2x2CellEta()       const {return f2x2CellEta ; }
-  Int_t    GetnxnCellEta()       const {return fnxnCellEta ; }
-  Int_t    Get2x2SuperModule()   const {return f2x2SM ; }
-  Int_t    GetnxnSuperModule()   const {return fnxnSM ; }
-
-  Int_t *  GetADCValuesLowGainMax2x2Sum()  {return fADCValuesLow2x2; }
-  Int_t *  GetADCValuesHighGainMax2x2Sum() {return fADCValuesHigh2x2; }
-  Int_t *  GetADCValuesLowGainMaxnxnSum()  {return fADCValuesLownxn; }
-  Int_t *  GetADCValuesHighGainMaxnxnSum() {return fADCValuesHighnxn; }
-
-  Float_t  GetL0Threshold() const           {return fL0Threshold ; } 
-  Float_t  GetL1JetLowPtThreshold() const   {return fL1JetLowPtThreshold ; }
-  Float_t  GetL1JetMediumPtThreshold()const {return fL1JetMediumPtThreshold ; }
-  Float_t  GetL1JetHighPtThreshold() const  {return fL1JetHighPtThreshold ; }
-
-  Float_t  GetPatchSize() const  {return fPatchSize ; }
-  Bool_t   IsSimulation() const {return fSimulation ; }
-  
+  Float_t  Get2x2MaxAmplitude()  const { return f2x2MaxAmp ; }
+  Float_t  GetnxnMaxAmplitude()  const { return fnxnMaxAmp ; }
+  Int_t    Get2x2CellPhi()       const { return f2x2CellPhi ; }
+  Int_t    GetnxnCellPhi()       const { return fnxnCellPhi ; }
+  Int_t    Get2x2CellEta()       const { return f2x2CellEta ; }
+  Int_t    GetnxnCellEta()       const { return fnxnCellEta ; }
+  Int_t    Get2x2SuperModule()   const { return f2x2SM ; }
+  Int_t    GetnxnSuperModule()   const { return fnxnSM ; }
+
+  Int_t *  GetADCValuesLowGainMax2x2Sum()  { return fADCValuesLow2x2; }
+  Int_t *  GetADCValuesHighGainMax2x2Sum() { return fADCValuesHigh2x2; }
+  Int_t *  GetADCValuesLowGainMaxnxnSum()  { return fADCValuesLownxn; }
+  Int_t *  GetADCValuesHighGainMaxnxnSum() { return fADCValuesHighnxn; }
+
+  Float_t  GetL0Threshold() const            { return fL0Threshold ; } 
+  Float_t  GetL1JetLowPtThreshold()    const { return fL1JetLowPtThreshold ; }
+  Float_t  GetL1JetMediumPtThreshold() const { return fL1JetMediumPtThreshold ; }
+  Float_t  GetL1JetHighPtThreshold()   const { return fL1JetHighPtThreshold ; }
+
+  Int_t    GetPatchSize()              const { return fPatchSize ; }
+  Int_t    GetIsolPatchSize()          const { return fIsolPatchSize ; }
+
+  Float_t  Get2x2AmpOutOfPatch()       const { return  f2x2AmpOutOfPatch ; }
+  Float_t  GetnxnAmpOutOfPatch()       const { return  fnxnAmpOutOfPatch ; }
+  Float_t  Get2x2AmpOutOfPatchThres()  const { return  f2x2AmpOutOfPatchThres ; }
+  Float_t  GetnxnAmpOutOfPatchThres()  const { return  fnxnAmpOutOfPatchThres ; } 
+
+  Bool_t   Is2x2Isol()                 const { return  fIs2x2Isol ; }
+  Bool_t   IsnxnIsol()                 const { return  fIsnxnIsol ; }
+
+  Bool_t   IsSimulation()              const { return fSimulation ; }
+  Bool_t   IsIsolatedInSuperModule()   const { return fIsolateInSuperModule ; }
+
   //Setters
   void     SetDigitsList(TClonesArray * digits)          
    {fDigitsList  = digits ; }
@@ -90,25 +108,31 @@ class AliEMCALTrigger : public AliTriggerDetector {
     {fL1JetHighPtThreshold   = amp; }
 
   void SetPatchSize(Int_t ps)                {fPatchSize = ps ; }
+  void SetIsolPatchSize(Int_t ps)          {fIsolPatchSize = ps ; }
+  void Set2x2AmpOutOfPatchThres(Float_t th) { f2x2AmpOutOfPatchThres = th; }
+  void SetnxnAmpOutOfPatchThres(Float_t th) { fnxnAmpOutOfPatchThres = th; }
   void SetSimulation(Bool_t sim )          {fSimulation = sim ; }
+  void SetIsolateInSuperModule(Bool_t isol )          {fIsolateInSuperModule = isol ; }
 
  private:
-  void MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus,const Int_t supermod, TMatrixD *ampmax2, TMatrixD *ampmaxn, AliEMCALGeometry * geom) ; 
+
+  Bool_t IsPatchIsolated(Int_t iPatchType, const TClonesArray * ampmods, const Int_t imod, const Int_t mtru, const Float_t maxamp, const Int_t maxphi, const Int_t maxeta) ;
+  
+  void MakeSlidingCell(const TClonesArray * amptrus, const TClonesArray * timeRtrus,const Int_t supermod, TMatrixD *ampmax2, TMatrixD *ampmaxn) ; 
+  
+  void SetTriggers(const TClonesArray * amptrus,const Int_t iSM, const TMatrixD *ampmax2, const TMatrixD *ampmaxn, const AliEMCALGeometry * geom) ;
   
 
-  void SetTriggers(const Int_t iSM, const TMatrixD *ampmax2, const TMatrixD *ampmaxn, AliEMCALGeometry *geom) ;
-    
  private: 
 
-  Float_t f2x2MaxAmp ;        //! Maximum 2x2 added amplitude (not overlapped) 
-  Int_t   f2x2CellPhi ;       //! upper right cell, row(phi)   
-  Int_t   f2x2CellEta ;       //! and column(eta)  
-  Int_t   f2x2SM ;            //! Super Module where maximum is found
-  Float_t fnxnMaxAmp ;        //! Maximum nxn added amplitude (overlapped)
-  Int_t   fnxnCellPhi ;       //! upper right cell, row(phi)   
-  Int_t   fnxnCellEta ;       //! and column(eta)
-  Int_t   fnxnSM ;            //! Super Module where maximum is found
+  Float_t f2x2MaxAmp ;         //! Maximum 2x2 added amplitude (not overlapped) 
+  Int_t   f2x2CellPhi ;        //! upper right cell, row(phi)   
+  Int_t   f2x2CellEta ;        //! and column(eta)  
+  Int_t   f2x2SM ;             //! Super Module where maximum is found
+  Float_t fnxnMaxAmp ;         //! Maximum nxn added amplitude (overlapped)
+  Int_t   fnxnCellPhi ;        //! upper right cell, row(phi)   
+  Int_t   fnxnCellEta ;        //! and column(eta)
+  Int_t   fnxnSM ;             //! Super Module where maximum is found
 
   Int_t*   fADCValuesHighnxn ; //! Sampled ADC high gain values for the nxn crystals amplitude sum
   Int_t*   fADCValuesLownxn  ; //! " low gain  " 
@@ -122,9 +146,28 @@ class AliEMCALTrigger : public AliTriggerDetector {
   Float_t fL1JetMediumPtThreshold ; //! L1 Medium pT trigger energy threshold
   Float_t fL1JetHighPtThreshold ;   //! L1 High pT trigger energy threshold
 
-  Int_t fPatchSize;                 //! Trigger patch factor, to be multiplied to 2x2 cells
-                                            // 0 means 2x2, 1 means nxn, 2 means 8x8 ...
+  Int_t   fNTRU;             //! Number of TRU per SuperModule (3)
+  Int_t   fNTRUEta ;         //! Number of crystal rows per Eta in one TRU (3)
+  Int_t   fNTRUPhi ;         //! Number of crystal rows per Phi in one TRU (1)
+  Int_t   fNCellsPhi;        //! Number of rows in a TRU (24)
+  Int_t   fNCellsEta;        //! Number of columns in a TRU (16)
+
+  Int_t fPatchSize;          //! Trigger patch factor, to be multiplied to 2x2 cells
+                             //  0 means 2x2, 1 means 4x4, 2 means 6x6 ...
+  Int_t fIsolPatchSize ;     //  Isolation patch size, number of rows or columns to add to 
+                             //  the 2x2 or nxn maximum amplitude patch. 
+                             //  1 means a patch around max amplitude of 2x2 of 4x4 and around         
+                             //  max ampl patch of 4x4 of 8x8 
+    
+  Float_t f2x2AmpOutOfPatch;      //  Amplitude in isolation cone minus maximum amplitude of the reference patch
+  Float_t fnxnAmpOutOfPatch; 
+  Float_t f2x2AmpOutOfPatchThres; //  Threshold to select a trigger as isolated on f2x2AmpOutOfPatch value
+  Float_t fnxnAmpOutOfPatchThres; 
+  Float_t fIs2x2Isol;             //  Patch is isolated if f2x2AmpOutOfPatchThres threshold is passed
+  Float_t fIsnxnIsol ; 
+
   Bool_t  fSimulation ;           //! Flag to do the trigger during simulation or reconstruction
+  Bool_t  fIsolateInSuperModule;  //! Flag to isolate trigger patch in SuperModule or in TRU acceptance
 
   ClassDef(AliEMCALTrigger,1)
 } ;