EMCAL/AliEMCALClusterizerv2.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 "AliEMCALClusterizerv2.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(AliEMCALClusterizerv2)
  55
  56 //____________________________________________________________________________
  57 AliEMCALClusterizerv2::AliEMCALClusterizerv2()
  58   : AliEMCALClusterizerv1(), fDoEnGradCut(1)
  59 {
  60   // ctor with the indication of the file where header Tree and digits Tree are stored
  61 }
  62
  63 //____________________________________________________________________________
  64 AliEMCALClusterizerv2::AliEMCALClusterizerv2(AliEMCALGeometry* geometry)
  65   : AliEMCALClusterizerv1(geometry), fDoEnGradCut(1)
  66 {
  67   // ctor with the indication of the file where header Tree and digits Tree are stored
  68   // use this contructor to avoid usage of Init() which uses runloader
  69   // change needed by HLT - MP
  70 }
  71
  72 //____________________________________________________________________________
  73 AliEMCALClusterizerv2::AliEMCALClusterizerv2(AliEMCALGeometry* geometry, AliEMCALCalibData* calib, AliCaloCalibPedestal* caloped)
  74   : AliEMCALClusterizerv1(geometry, calib, caloped), fDoEnGradCut(1)
  75 {
  76   // ctor, geometry and calibration are initialized elsewhere.
  77 }
  78
  79 //____________________________________________________________________________
  80 AliEMCALClusterizerv2::~AliEMCALClusterizerv2()
  81 {
  82   // dtor
  83 }
  84
  85 //____________________________________________________________________________
  86 Int_t AliEMCALClusterizerv2::AreNeighbours(AliEMCALDigit* d1, AliEMCALDigit* d2, Bool_t& shared) const
  87 {
  88   // Gives the neighbourness of two digits = 0 are not neighbour; continue searching
  89   //                                       = 1 are neighbour
  90   //                                       = 2 is in different SM; continue searching
  91   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
  92   // neighbours are defined as digits having at least a common side
  93   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
  94   //                                      which is compared to a digit (d2)  not yet in a cluster
  95
  96   if (fDoEnGradCut) {
  97     if (d2->GetCalibAmp()>d1->GetCalibAmp())
  98       return 3; // energy of neighboring cell should be smaller in order to become a neighbor
  99   }
 100   return AliEMCALClusterizerv1::AreNeighbours(d1,d2,shared);
 101 }
 102
 103 //____________________________________________________________________________
 104 void AliEMCALClusterizerv2::MakeClusters()
 105 {
 106   // Make list of clusters. Start from highest energy cell.
 107
 108   if (fGeom==0)
 109     AliFatal("Did not get geometry from EMCALLoader");
 110
 111   fRecPoints->Delete();
 112
 113   // set up TObjArray with pointers to digits to work on, calibrate digits
 114   TObjArray digitsC;
 115   Double_t ehs = 0.0;
 116   AliEMCALDigit *digit = 0;
 117   TIter nextdigit(fDigitsArr);
 118   while ( (digit = static_cast<AliEMCALDigit*>(nextdigit())) ) {
 119     Float_t dEnergyCalibrated = digit->GetAmplitude();
 120     Float_t time              = digit->GetTime();
 121     Calibrate(dEnergyCalibrated, time ,digit->GetId());
 122     digit->SetCalibAmp(dEnergyCalibrated);
 123     digit->SetTime(time);
 124     if (dEnergyCalibrated < fMinECut)
 125       continue;
 126     if (!fGeom->CheckAbsCellId(digit->GetId()))
 127       continue;
 128     ehs += dEnergyCalibrated;
 129     digitsC.AddLast(digit);
 130   }
 131
 132   AliDebug(1,Form("MakeClusters: Number of digits %d  -> ehs %f (minE %f)\n",
 133                   fDigitsArr->GetEntries(),ehs,fMinECut));
 134
 135   TIter nextdigitC(&digitsC);
 136   while (1) {
 137     Int_t   iMaxEnergyDigit = -1;
 138     Float_t dMaxEnergyDigit = -1;
 139     AliEMCALDigit *pMaxEnergyDigit = 0;
 140     nextdigitC.Reset();
 141     while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) ) {
 142       Float_t dEnergyCalibrated = digit->GetCalibAmp();
 143       if (dEnergyCalibrated>fECAClusteringThreshold && dEnergyCalibrated>dMaxEnergyDigit) {
 144         dMaxEnergyDigit = dEnergyCalibrated;
 145         iMaxEnergyDigit = digit->GetId();
 146         pMaxEnergyDigit = digit;
 147       }
 148     }
 149     if (iMaxEnergyDigit<0 || digitsC.GetEntries() <= 0) {
 150       break;
 151     }
 152
 153     if (fNumberOfECAClusters>=fRecPoints->GetSize())
 154       fRecPoints->Expand(2*fNumberOfECAClusters+1);
 155
 156     AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("");
 157     recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 158     recPoint->AddDigit(*pMaxEnergyDigit, pMaxEnergyDigit->GetCalibAmp(), kFALSE);
 159     fRecPoints->AddAt(recPoint, fNumberOfECAClusters++);
 160     digitsC.Remove(pMaxEnergyDigit);
 161     TObjArray clusterDigits;
 162     clusterDigits.AddLast(pMaxEnergyDigit);
 163     TIter nextClusterDigit(&clusterDigits);
 164     Float_t time = pMaxEnergyDigit->GetTime();
 165
 166     AliDebug(1,Form("MakeClusters: Max digit found id = %d, ene = %f , clus.th. = %f \n",
 167                     iMaxEnergyDigit, dMaxEnergyDigit, fECAClusteringThreshold));
 168
 169     while ( (digit = static_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
 170       TIter nextdigitN(&digitsC);
 171       AliEMCALDigit *digitN = 0; // digi neighbor
 172       while ( (digitN = static_cast<AliEMCALDigit*>(nextdigitN())) ) { // scan over all digits to look for neighbours
 173         //Do not add digits with too different time
 174         if (TMath::Abs(time - digitN->GetTime()) > fTimeCut )
 175           continue;
 176         Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
 177         if (AreNeighbours(digit, digitN, shared)==1) {
 178           recPoint->AddDigit(*digitN, digitN->GetCalibAmp(), shared);
 179           clusterDigits.AddLast(digitN);
 180           digitsC.Remove(digitN);
 181         }
 182       }
 183     }
 184     AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
 185   }
 186   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
 187 }
 188