EMCAL/AliEMCALClusterizerNxN.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 // This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
  17 // Algorithm:
  18 // 1. peek the most energetic cell
  19 // 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
  20 // 3. remove the cells contributing to the cluster
  21 // 4. start from 1 for the remaining clusters
  22 // 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
  23 // - for high energy clusters check the surrounding of the 3x3 clusters for extra energy
  24 // (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
  25 // Use Case:
  26 //  root [0] AliEMCALClusterizerNxN * cl = new AliEMCALClusterizerNxN("galice.root")
  27 //  Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
  28 //               //reads gAlice from header file "..."
  29 //  root [1] cl->ExecuteTask()
  30 //               //finds RecPoints in all events stored in galice.root
  31 //  root [2] cl->SetDigitsBranch("digits2")
  32 //               //sets another title for Digitis (input) branch
  33 //  root [3] cl->SetRecPointsBranch("recp2")
  34 //               //sets another title four output branches
  35 //  root [4] cl->SetTowerLocalMaxCut(0.03)
  36 //               //set clusterization parameters
  37 //  root [5] cl->ExecuteTask("deb all time")
  38 //               //once more finds RecPoints options are
  39 //               // deb - print number of found rec points
  40 //               // deb all - print number of found RecPoints and some their characteristics
  41 //               // time - print benchmarking results
  42
  43 // --- ROOT system ---
  44 #include <TMath.h>
  45 #include <TMinuit.h>
  46 #include <TTree.h>
  47 #include <TBenchmark.h>
  48 #include <TBrowser.h>
  49 #include <TROOT.h>
  50 #include <TClonesArray.h>
  51
  52 // --- Standard library ---
  53 #include <cassert>
  54
  55 // --- AliRoot header files ---
  56 #include "AliLog.h"
  57 #include "AliEMCALClusterizerNxN.h"
  58 #include "AliEMCALRecPoint.h"
  59 #include "AliEMCALDigit.h"
  60 #include "AliEMCALGeometry.h"
  61 #include "AliCaloCalibPedestal.h"
  62 #include "AliEMCALCalibData.h"
  63 #include "AliESDCaloCluster.h"
  64 #include "AliEMCALUnfolding.h"
  65
  66 ClassImp(AliEMCALClusterizerNxN)
  67
  68 //____________________________________________________________________________
  69 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
  70   : AliEMCALClusterizer(), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
  71 {
  72   // ctor with the indication of the file where header Tree and digits Tree are stored
  73 }
  74
  75 //____________________________________________________________________________
  76 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
  77   : AliEMCALClusterizer(geometry), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
  78 {
  79   // ctor with the indication of the file where header Tree and digits Tree are stored
  80   // use this contructor to avoid usage of Init() which uses runloader
  81   // change needed by HLT - MP
  82 }
  83
  84 //____________________________________________________________________________
  85 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
  86 : AliEMCALClusterizer(geometry, calib, caloped), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0)
  87 {
  88   // ctor, geometry and calibration are initialized elsewhere.
  89 }
  90
  91 //____________________________________________________________________________
  92 AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
  93 {
  94   // dtor
  95 }
  96
  97 //____________________________________________________________________________
  98 void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
  99 {
 100   // Steering method to perform clusterization for the current event
 101   // in AliEMCALLoader
 102
 103   if (strstr(option,"tim"))
 104     gBenchmark->Start("EMCALClusterizer");
 105
 106   if (strstr(option,"print"))
 107     Print("");
 108
 109   //Get calibration parameters from file or digitizer default values.
 110   GetCalibrationParameters();
 111
 112   //Get dead channel map from file or digitizer default values.
 113   GetCaloCalibPedestal();
 114
 115   MakeClusters();  //only the real clusters
 116
 117   if (fToUnfold) {
 118     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
 119     fClusterUnfolding->MakeUnfolding();
 120   }
 121
 122   //Evaluate position, dispersion and other RecPoint properties for EC section
 123   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
 124     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 125     if (rp) {
 126       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
 127       AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
 128       //For each rec.point set the distance to the nearest bad crystal
 129       if (fCaloPed)
 130         rp->EvalDistanceToBadChannels(fCaloPed);
 131     }
 132   }
 133
 134   fRecPoints->Sort();
 135
 136   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
 137     AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 138     if (rp) {
 139       rp->SetIndexInList(index);
 140     }
 141     else AliFatal("RecPoint NULL!!");
 142   }
 143
 144   if (fTreeR)
 145     fTreeR->Fill();
 146
 147   if (strstr(option,"deb") || strstr(option,"all"))
 148     PrintRecPoints(option);
 149
 150   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
 151
 152   if (strstr(option,"tim")) {
 153     gBenchmark->Stop("EMCALClusterizer");
 154     printf("Exec took %f seconds for Clusterizing",
 155            gBenchmark->GetCpuTime("EMCALClusterizer"));
 156   }
 157 }
 158
 159 //____________________________________________________________________________
 160 Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
 161 {
 162   // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
 163   //                                       = 1 are neighbour
 164   //                                       = 2 is in different SM; continue searching
 165   // In case it is in different SM, but same phi rack, check if neigbours at eta=0
 166   // neighbours are defined as digits having at least a common side
 167   // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 168   //                                      which is compared to a digit (d2)  not yet in a cluster
 169
 170   if (fEnergyGrad) { //false by default
 171     if (d2->GetCalibAmp()>d1->GetCalibAmp())
 172       return 3; // energy of neighboring cell should be smaller in order to become a neighbor
 173   }
 174
 175   Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 176   Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 177   Int_t rowdiff=0, coldiff=0;
 178
 179   shared = kFALSE;
 180
 181   fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
 182   fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
 183   fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
 184   fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
 185
 186   //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
 187   if (nSupMod1 != nSupMod2 )
 188     {
 189       //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
 190       Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
 191       Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
 192
 193       if (!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
 194
 195       // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
 196       // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 197       if (nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
 198       else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
 199
 200       shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
 201
 202     }//Different SM, same phi
 203
 204   rowdiff = TMath::Abs(iphi1 - iphi2);
 205   coldiff = TMath::Abs(ieta1 - ieta2);
 206
 207   // neighbours +-1 in col and row
 208   if ( TMath::Abs(coldiff) <= fNColDiff && TMath::Abs(rowdiff) <= fNRowDiff)
 209     {
 210
 211       AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
 212                        d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
 213
 214       return 1;
 215     }//Neighbours
 216   else
 217     {
 218       AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
 219                        d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
 220       shared = kFALSE;
 221       return 2;
 222     }//Not neighbours
 223 }
 224
 225 //____________________________________________________________________________
 226 void AliEMCALClusterizerNxN::MakeClusters()
 227 {
 228   // Make clusters
 229
 230   if (fGeom==0)
 231     AliFatal("Did not get geometry from EMCALLoader");
 232
 233   fNumberOfECAClusters = 0;
 234   fRecPoints->Delete();
 235
 236   // Set up TObjArray with pointers to digits to work on, calibrate digits
 237   TObjArray digitsC;
 238   TIter nextdigit(fDigitsArr);
 239   AliEMCALDigit *digit = 0;
 240   while ( (digit = static_cast<AliEMCALDigit*>(nextdigit())) ) {
 241     Float_t dEnergyCalibrated = digit->GetAmplitude();
 242     Float_t time              = digit->GetTime();
 243     Calibrate(dEnergyCalibrated,time ,digit->GetId());
 244     digit->SetCalibAmp(dEnergyCalibrated);
 245     digit->SetTime(time);
 246     digitsC.AddLast(digit);
 247   }
 248
 249   TIter nextdigitC(&digitsC);
 250
 251   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f)\n",
 252                   fDigitsArr->GetEntries(),fMinECut));
 253
 254   Bool_t bDone = kFALSE;
 255   while ( bDone != kTRUE )
 256   {
 257     //first sort the digits:
 258     Int_t iMaxEnergyDigit = -1;
 259     Float_t dMaxEnergyDigit = -1;
 260     AliEMCALDigit *pMaxEnergyDigit = 0;
 261     nextdigitC.Reset();
 262     while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) )
 263     { // scan over the list of digitsC
 264       Float_t dEnergyCalibrated = digit->GetCalibAmp();
 265
 266       if (fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fMinECut) // no threshold by default!
 267       {                                                                          // needs to be set in OCDB!
 268         if (dEnergyCalibrated > dMaxEnergyDigit)
 269         {
 270           dMaxEnergyDigit = dEnergyCalibrated;
 271           iMaxEnergyDigit = digit->GetId();
 272           pMaxEnergyDigit = digit;
 273         }
 274       }
 275     }
 276     if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0)
 277     {
 278       bDone = kTRUE;
 279       continue;
 280     }
 281
 282     AliDebug (2, Form("Max digit found: %1.5f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
 283
 284     // keep the candidate digits in a list
 285     TList clusterDigitList;
 286     clusterDigitList.SetOwner(kFALSE);
 287     clusterDigitList.AddLast(pMaxEnergyDigit);
 288
 289     Double_t clusterCandidateEnergy = dMaxEnergyDigit;
 290
 291     // now loop over the rest of the digits and cluster into NxN cluster
 292     // we do not actually cluster yet: we keep them in the list clusterDigitList
 293     nextdigitC.Reset();
 294     while ( (digit = static_cast<AliEMCALDigit *>(nextdigitC())) )
 295     { // scan over the list of digitsC
 296       if (digit == pMaxEnergyDigit) continue;
 297       Float_t dEnergyCalibrated = digit->GetCalibAmp();
 298       AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
 299       if (fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0  )
 300       {
 301         Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
 302         if (TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
 303         Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
 304         if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!!
 305         {
 306           clusterDigitList.AddLast(digit);
 307           clusterCandidateEnergy += dEnergyCalibrated;
 308         }
 309       }
 310     }// loop over the next digits
 311
 312     // start a cluster here only if a cluster energy is larger than clustering threshold
 313     AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
 314     if (clusterCandidateEnergy > fECAClusteringThreshold)
 315     {
 316       if (fNumberOfECAClusters >= fRecPoints->GetSize())
 317         fRecPoints->Expand(2*fNumberOfECAClusters+1);
 318
 319       (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ;
 320       AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(fNumberOfECAClusters) ) ;
 321       // AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("");
 322       // fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
 323       // recPoint = static_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters));
 324       if (recPoint) {
 325         fNumberOfECAClusters++;
 326         recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 327
 328         AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
 329         for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
 330         {
 331           digit = (AliEMCALDigit*)clusterDigitList.At(idig);
 332           AliDebug(5, Form(" Adding digit %d", digit->GetId()));
 333           // note: this way the sharing info is lost!
 334           recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
 335           digitsC.Remove(digit);
 336         }
 337       }// recpoint
 338     }
 339     else
 340     {
 341       // we do not want to start clustering in the same spot!
 342       // but in this case we may NOT reuse this seed for another cluster!
 343       // need a better bookeeping?
 344       digitsC.Remove(pMaxEnergyDigit);
 345     }
 346
 347     AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));
 348   } // while ! done
 349
 350   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
 351 }