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 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
  75 : AliEMCALClusterizer(geometry, calib, caloped)
  76 {
  77   // ctor, geometry and calibration are initialized elsewhere.
  78 }
  79
  80 //____________________________________________________________________________
  81   AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
  82 {
  83   // dtor
  84 }
  85
  86 //____________________________________________________________________________
  87 void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
  88 {
  89   // Steering method to perform clusterization for the current event
  90   // in AliEMCALLoader
  91
  92   if(strstr(option,"tim"))
  93     gBenchmark->Start("EMCALClusterizer");
  94
  95   if(strstr(option,"print"))
  96     Print("");
  97
  98   //Get calibration parameters from file or digitizer default values.
  99   GetCalibrationParameters();
 100
 101   //Get dead channel map from file or digitizer default values.
 102   GetCaloCalibPedestal();
 103
 104   fNumberOfECAClusters = 0;
 105
 106   MakeClusters();  //only the real clusters
 107
 108   if(fToUnfold){
 109     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
 110     fClusterUnfolding->MakeUnfolding();
 111   }
 112
 113   //Evaluate position, dispersion and other RecPoint properties for EC section
 114   Int_t index;
 115   for(index = 0; index < fRecPoints->GetEntries(); index++) {
 116     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 117     if(rp){
 118       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
 119       //For each rec.point set the distance to the nearest bad crystal
 120       if (fCaloPed)
 121         rp->EvalDistanceToBadChannels(fCaloPed);
 122     }
 123     else AliFatal("Null rec point in list!");
 124   }
 125
 126   fRecPoints->Sort();
 127
 128   for(index = 0; index < fRecPoints->GetEntries(); index++) {
 129     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 130     if(rp){
 131       rp->SetIndexInList(index);
 132       rp->Print();
 133     }
 134     else AliFatal("Null rec point in list!");
 135   }
 136
 137   if (fTreeR)
 138     fTreeR->Fill();
 139
 140   if(strstr(option,"deb") || strstr(option,"all"))
 141     PrintRecPoints(option);
 142
 143   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
 144
 145   if(strstr(option,"tim")){
 146     gBenchmark->Stop("EMCALClusterizer");
 147     printf("Exec took %f seconds for Clusterizing",
 148            gBenchmark->GetCpuTime("EMCALClusterizer"));
 149   }
 150 }
 151
 152 //____________________________________________________________________________
 153 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
 154 {
 155   // Gives the neighbourness of two digits = 0 are not neighbour; continue searching
 156   //                                       = 1 are neighbour
 157   //                                       = 2 is in different SM; continue searching
 158   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
 159   // neighbours are defined as digits having at least a common side
 160   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 161   //                                      which is compared to a digit (d2)  not yet in a cluster
 162
 163   Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 164   Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 165
 166   shared = kFALSE;
 167
 168   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
 169   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
 170   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
 171   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
 172
 173   //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
 174   if (nSupMod1 != nSupMod2 ) {
 175     //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
 176     Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
 177     Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
 178
 179     if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
 180
 181     // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
 182     // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 183     if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
 184     else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
 185
 186     shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
 187   } //Different SM, same phi
 188
 189   Int_t rowdiff = TMath::Abs(iphi1 - iphi2);
 190   Int_t coldiff = TMath::Abs(ieta1 - ieta2);
 191
 192   // neighbours with at least common side; May 11, 2007
 193   if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {
 194     //Diagonal?
 195     //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
 196
 197     if (gDebug == 2)
 198       printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
 199              d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);
 200     return 1;
 201   } //Neighbours
 202   else {
 203     shared = kFALSE;
 204     return 2;
 205   } //Not neighbours
 206 }
 207
 208 //____________________________________________________________________________
 209 void AliEMCALClusterizerv1::MakeClusters()
 210 {
 211   // Steering method to construct the clusters stored in a list of Reconstructed Points
 212   // A cluster is defined as a list of neighbour digits
 213   // Mar 03, 2007 by PAI
 214
 215   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
 216
 217   fRecPoints->Delete();
 218
 219   // Set up TObjArray with pointers to digits to work on calibrated digits
 220   TObjArray *digitsC = new TObjArray();
 221   AliEMCALDigit *digit;
 222   Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
 223   TIter nextdigit(fDigitsArr);
 224   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // calibrate and clean up digits
 225     dEnergyCalibrated =  digit->GetAmplitude();
 226     time              =  digit->GetTime();
 227     Calibrate(dEnergyCalibrated,time,digit->GetId());
 228     digit->SetCalibAmp(dEnergyCalibrated);
 229     digit->SetTime(time);
 230     if ( dEnergyCalibrated < fMinECut){
 231       continue;
 232     }
 233     else if (!fGeom->CheckAbsCellId(digit->GetId()))
 234       continue;
 235     else{
 236       ehs += dEnergyCalibrated;
 237       digitsC->AddLast(digit);
 238     }
 239   }
 240
 241   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
 242                   fDigitsArr->GetEntries(),fMinECut,ehs));
 243
 244   TIter nextdigitC(digitsC);
 245   while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
 246     TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
 247     dEnergyCalibrated = digit->GetCalibAmp();
 248     time              = digit->GetTime();
 249     if(fGeom->CheckAbsCellId(digit->GetId()) && ( dEnergyCalibrated > 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         recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
 261         TObjArray clusterDigits;
 262         clusterDigits.AddLast(digit);
 263         digitsC->Remove(digit);
 264
 265         AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(), dEnergyCalibrated, fECAClusteringThreshold));  //Time or TimeR?
 266
 267         // Grow cluster by finding neighbours
 268         TIter nextClusterDigit(&clusterDigits);
 269
 270         while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
 271           TIter nextdigitN(digitsC);
 272           AliEMCALDigit *digitN = 0; // digi neighbor
 273           while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
 274             //Do not add digits with too different time
 275             Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
 276             if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
 277             if (AreNeighbours(digit, digitN, shared)==1) {      // call (digit,digitN) in THAT order !!!!!
 278               recPoint->AddDigit(*digitN, digitN->GetCalibAmp(), shared);
 279               clusterDigits.AddLast(digitN);
 280               digitsC->Remove(digitN);
 281             } // if(ineb==1)
 282           } // scan over digits
 283         } // scan over digits already in cluster
 284
 285         AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
 286       }//recpoint
 287       else AliFatal("Null recpoint in array!");
 288     } // If seed found
 289   } // while digit
 290
 291   delete digitsC;
 292
 293   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
 294 }