EMCAL/AliEMCALClusterizerv1.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
  16 /* $Id$ */
  17
  18 //-- Author: Yves Schutz (SUBATECH)  & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
  19 //--         Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
  20 //////////////////////////////////////////////////////////////////////////////
  21 //  Clusterization class. Performs clusterization (collects neighbouring active cells) and
  22 //  unfolds the clusters having several local maxima.
  23 //  Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
  24 //  EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
  25 //  parameters including input digits branch title, thresholds etc.)
  26 //
  27
  28 // --- ROOT system ---
  29
  30 #include <TFile.h>
  31 #include <TMath.h>
  32 #include <TMinuit.h>
  33 #include <TTree.h>
  34 #include <TBenchmark.h>
  35 #include <TBrowser.h>
  36 #include <TROOT.h>
  37 #include <TList.h>
  38 #include <TClonesArray.h>
  39
  40 // --- Standard library ---
  41 #include <cassert>
  42
  43 // --- AliRoot header files ---
  44 #include "AliLog.h"
  45 #include "AliEMCALClusterizerv1.h"
  46 #include "AliEMCALRecPoint.h"
  47 #include "AliEMCALDigit.h"
  48 #include "AliEMCALGeometry.h"
  49 #include "AliCaloCalibPedestal.h"
  50 #include "AliEMCALCalibData.h"
  51 #include "AliESDCaloCluster.h"
  52 #include "AliEMCALUnfolding.h"
  53
  54 ClassImp(AliEMCALClusterizerv1)
  55
  56 //____________________________________________________________________________
  57 AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer()
  58 {
  59   // ctor with the indication of the file where header Tree and digits Tree are stored
  60
  61   Init() ;
  62 }
  63
  64 //____________________________________________________________________________
  65 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
  66   : AliEMCALClusterizer(geometry)
  67 {
  68   // ctor with the indication of the file where header Tree and digits Tree are stored
  69   // use this contructor to avoid usage of Init() which uses runloader
  70   // change needed by HLT - MP
  71
  72 }
  73
  74 //____________________________________________________________________________
  75 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
  76 : AliEMCALClusterizer(geometry, calib, caloped)
  77 {
  78         // ctor, geometry and calibration are initialized elsewhere.
  79
  80 }
  81
  82
  83 //____________________________________________________________________________
  84   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
  85 {
  86   // dtor
  87 }
  88
  89 //____________________________________________________________________________
  90 void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
  91 {
  92   // Steering method to perform clusterization for the current event
  93   // in AliEMCALLoader
  94
  95   if(strstr(option,"tim"))
  96     gBenchmark->Start("EMCALClusterizer");
  97
  98   if(strstr(option,"print"))
  99     Print("") ;
 100
 101   //Get calibration parameters from file or digitizer default values.
 102   GetCalibrationParameters() ;
 103
 104   //Get dead channel map from file or digitizer default values.
 105   GetCaloCalibPedestal() ;
 106
 107   fNumberOfECAClusters = 0;
 108
 109   MakeClusters() ;  //only the real clusters
 110
 111   if(fToUnfold){
 112     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
 113     fClusterUnfolding->MakeUnfolding();
 114   }
 115
 116   //Evaluate position, dispersion and other RecPoint properties for EC section
 117   Int_t index ;
 118   for(index = 0; index < fRecPoints->GetEntries(); index++) {
 119     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 120     if(rp){
 121       rp->EvalAll(fECAW0,fDigitsArr) ;
 122       //For each rec.point set the distance to the nearest bad crystal
 123             rp->EvalDistanceToBadChannels(fCaloPed);
 124     }
 125     else AliFatal("Null rec point in list!");
 126   }
 127
 128   fRecPoints->Sort() ;
 129
 130   for(index = 0; index < fRecPoints->GetEntries(); index++) {
 131     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 132     if(rp){
 133       rp->SetIndexInList(index) ;
 134       rp->Print();
 135     }
 136     else AliFatal("Null rec point in list!");
 137   }
 138
 139   fTreeR->Fill();
 140
 141   if(strstr(option,"deb") || strstr(option,"all"))
 142     PrintRecPoints(option) ;
 143
 144   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
 145
 146   fRecPoints->Delete();
 147
 148   if(strstr(option,"tim")){
 149     gBenchmark->Stop("EMCALClusterizer");
 150     printf("Exec took %f seconds for Clusterizing",
 151            gBenchmark->GetCpuTime("EMCALClusterizer"));
 152   }
 153 }
 154
 155 //____________________________________________________________________________
 156 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
 157 {
 158   // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
 159   //                                       = 1 are neighbour
 160   //                                       = 2 is in different SM; continue searching
 161   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
 162   // neighbours are defined as digits having at least a common side
 163   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 164   //                                      which is compared to a digit (d2)  not yet in a cluster
 165
 166   static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 167   static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 168
 169   shared = kFALSE;
 170
 171   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
 172   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
 173   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
 174   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
 175
 176   //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
 177   if(nSupMod1 != nSupMod2 ) {
 178     //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
 179     Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
 180     Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
 181
 182     if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
 183
 184     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
 185     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 186     if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
 187     else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
 188
 189     shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
 190
 191   }//Different SM, same phi
 192
 193   Int_t rowdiff = TMath::Abs(iphi1 - iphi2);
 194   Int_t coldiff = TMath::Abs(ieta1 - ieta2) ;
 195
 196   // neighbours with at least common side; May 11, 2007
 197   if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {
 198     //Diagonal?
 199     //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
 200
 201     if (gDebug == 2)
 202       printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
 203              d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);
 204
 205     return 1;
 206   }//Neighbours
 207   else {
 208     shared = kFALSE;
 209     return 2 ;
 210   }//Not neighbours
 211 }
 212
 213 //____________________________________________________________________________
 214 void AliEMCALClusterizerv1::MakeClusters()
 215 {
 216   // Steering method to construct the clusters stored in a list of Reconstructed Points
 217   // A cluster is defined as a list of neighbour digits
 218   // Mar 03, 2007 by PAI
 219
 220   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
 221
 222   fRecPoints->Clear();
 223
 224   // Set up TObjArray with pointers to digits to work on
 225   TObjArray *digitsC = new TObjArray();
 226   TIter nextdigit(fDigitsArr);
 227   AliEMCALDigit *digit;
 228   while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
 229     digitsC->AddLast(digit);
 230   }
 231
 232   double e = 0.0, ehs = 0.0;
 233   TIter nextdigitC(digitsC);
 234   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
 235     e = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());//Time or TimeR?
 236     if ( e < fMinECut) //|| digit->GetTimeR() > fTimeCut ) // time window of cell checked in calibrate
 237       digitsC->Remove(digit);
 238     else
 239       ehs += e;
 240   }
 241   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
 242                   fDigitsArr->GetEntries(),fMinECut,ehs));
 243
 244   nextdigitC.Reset();
 245
 246   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
 247     TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
 248
 249     if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()) > fECAClusteringThreshold  ) ){
 250       // start a new Tower RecPoint
 251       if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
 252
 253       AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ;
 254       fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
 255       recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
 256       if(recPoint){
 257         fNumberOfECAClusters++ ;
 258
 259         recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 260
 261         recPoint->AddDigit(*digit, Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()),kFALSE) ; //Time or TimeR?
 262         TObjArray clusterDigits;
 263         clusterDigits.AddLast(digit);
 264         digitsC->Remove(digit) ;
 265
 266         AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(),
 267                         Calibrate(digit->GetAmplitude(),digit->GetTime(),digit->GetId()), fECAClusteringThreshold));  //Time or TimeR?
 268         Float_t time = digit->GetTime();//Time or TimeR?
 269         // Grow cluster by finding neighbours
 270         TIter nextClusterDigit(&clusterDigits);
 271         while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
 272           TIter nextdigitN(digitsC);
 273           AliEMCALDigit *digitN = 0; // digi neighbor
 274           while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
 275
 276             //Do not add digits with too different time
 277             Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
 278             if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
 279             if (AreNeighbours(digit, digitN, shared)==1) {      // call (digit,digitN) in THAT order !!!!!
 280               recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmplitude(), digitN->GetTime(), digitN->GetId()),shared) ;//Time or TimeR?
 281               clusterDigits.AddLast(digitN) ;
 282               digitsC->Remove(digitN) ;
 283             } // if(ineb==1)
 284           } // scan over digits
 285         } // scan over digits already in cluster
 286
 287         AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
 288       }//recpoint
 289       else AliFatal("Null recpoint in array!");
 290     } // If seed found
 291   } // while digit
 292
 293   delete digitsC ;
 294
 295   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
 296 }
 297