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
  65 ClassImp(AliEMCALClusterizerNxN)
  66
  67 Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE;
  68
  69 //____________________________________________________________________________
  70 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN()
  71   : AliEMCALClusterizer()
  72 {
  73   // ctor with the indication of the file where header Tree and digits Tree are stored
  74 }
  75
  76 //____________________________________________________________________________
  77 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry)
  78   : AliEMCALClusterizer(geometry)
  79 {
  80   // ctor with the indication of the file where header Tree and digits Tree are stored
  81   // use this contructor to avoid usage of Init() which uses runloader
  82   // change needed by HLT - MP
  83
  84 }
  85
  86 //____________________________________________________________________________
  87 AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
  88 : AliEMCALClusterizer(geometry, calib, caloped)
  89 {
  90         // ctor, geometry and calibration are initialized elsewhere.
  91
  92 }
  93
  94
  95 //____________________________________________________________________________
  96   AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN()
  97 {
  98   // dtor
  99 }
 100
 101
 102 //____________________________________________________________________________
 103 void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option)
 104 {
 105   // Steering method to perform clusterization for the current event
 106   // in AliEMCALLoader
 107
 108   if(strstr(option,"tim"))
 109     gBenchmark->Start("EMCALClusterizer");
 110
 111   if(strstr(option,"print"))
 112     Print("") ;
 113
 114   //Get calibration parameters from file or digitizer default values.
 115   GetCalibrationParameters() ;
 116
 117   //Get dead channel map from file or digitizer default values.
 118   GetCaloCalibPedestal() ;
 119
 120   fNumberOfECAClusters = 0;
 121
 122   MakeClusters() ;  //only the real clusters
 123
 124   if(fToUnfold)
 125     MakeUnfolding() ;
 126
 127   Int_t index ;
 128
 129   //Evaluate position, dispersion and other RecPoint properties for EC section
 130   for(index = 0; index < fRecPoints->GetEntries(); index++)
 131     {
 132       dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ;
 133       AliDebug(5, Form("MAX INDEX %d ", dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->GetMaximalEnergyIndex()));
 134       //For each rec.point set the distance to the nearest bad crystal
 135       dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index))->EvalDistanceToBadChannels(fCaloPed);
 136     }
 137
 138   fRecPoints->Sort() ;
 139
 140   for(index = 0; index < fRecPoints->GetEntries(); index++)
 141     {
 142       (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->SetIndexInList(index) ;
 143       (dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index)))->Print();
 144     }
 145
 146   fTreeR->Fill();
 147
 148   if(strstr(option,"deb") || strstr(option,"all"))
 149     PrintRecPoints(option) ;
 150
 151   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
 152
 153   fRecPoints->Delete();
 154
 155   if(strstr(option,"tim")){
 156     gBenchmark->Stop("EMCALClusterizer");
 157     printf("Exec took %f seconds for Clusterizing",
 158            gBenchmark->GetCpuTime("EMCALClusterizer"));
 159   }
 160 }
 161
 162 //____________________________________________________________________________
 163 Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
 164 {
 165         // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
 166         //                                       = 1 are neighbour
 167         //                                       = 2 is in different SM; continue searching
 168         // In case it is in different SM, but same phi rack, check if neigbours at eta=0
 169         // neighbours are defined as digits having at least a common side
 170         // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
 171         //                                      which is compared to a digit (d2)  not yet in a cluster
 172
 173         static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
 174         static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
 175         static Int_t rowdiff=0, coldiff=0;
 176
 177         shared = kFALSE;
 178
 179         fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
 180         fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
 181         fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
 182         fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
 183
 184         //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
 185         if(nSupMod1 != nSupMod2 )
 186           {
 187             //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
 188             Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
 189             Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
 190
 191             if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
 192
 193             // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
 194             // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
 195             if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
 196             else           ieta2+=AliEMCALGeoParams::fgkEMCALCols;
 197
 198             shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
 199
 200         }//Different SM, same phi
 201
 202         rowdiff = TMath::Abs(iphi1 - iphi2);
 203         coldiff = TMath::Abs(ieta1 - ieta2) ;
 204
 205         // neighbours +-1 in col and row
 206         if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2)
 207           {
 208
 209             AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
 210                              d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
 211
 212             return 1;
 213           }//Neighbours
 214         else
 215           {
 216             AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
 217                              d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared));
 218             shared = kFALSE;
 219             return 2 ;
 220           }//Not neighbours
 221 }
 222
 223 //____________________________________________________________________________
 224 void AliEMCALClusterizerNxN::MakeClusters()
 225 {
 226   // Steering method to construct the clusters stored in a list of Reconstructed Points
 227   // A cluster is defined as a list of neighbour digits
 228   // Mar 03, 2007 by PAI
 229
 230   if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
 231
 232   fRecPoints->Clear();
 233
 234   // Set up TObjArray with pointers to digits to work on
 235   //TObjArray *digitsC = new TObjArray();
 236   TObjArray digitsC;
 237   TIter nextdigit(fDigitsArr);
 238   AliEMCALDigit *digit = 0;
 239   while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
 240     digitsC.AddLast(digit);
 241   }
 242
 243   TIter nextdigitC(&digitsC);
 244
 245   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f)\n",
 246                   fDigitsArr->GetEntries(),fMinECut));
 247
 248   Bool_t bDone = kFALSE;
 249   while ( bDone != kTRUE )
 250     {
 251       //first sort the digits:
 252       Int_t iMaxEnergyDigit = -1;
 253       Float_t dMaxEnergyDigit = -1;
 254       AliEMCALDigit *pMaxEnergyDigit = 0;
 255       nextdigitC.Reset();
 256       while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
 257         { // scan over the list of digitsC
 258           Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
 259           //AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
 260
 261           //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold  )
 262           if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold!
 263             {
 264               if (dEnergyCalibrated > dMaxEnergyDigit)
 265                 {
 266                   dMaxEnergyDigit = dEnergyCalibrated;
 267                   iMaxEnergyDigit = digit->GetId();
 268                   pMaxEnergyDigit = digit;
 269                 }
 270             }
 271         }
 272
 273       if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0)
 274         {
 275           bDone = kTRUE;
 276           continue;
 277         }
 278
 279       AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit));
 280       AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit));
 281
 282       // keep the candidate digits in a list
 283       TList clusterDigitList;
 284       clusterDigitList.SetOwner(kFALSE);
 285       clusterDigitList.AddLast(pMaxEnergyDigit);
 286
 287       Double_t clusterCandidateEnergy = dMaxEnergyDigit;
 288
 289       // now loop over the resto of the digits and cluster into NxN cluster
 290       // we do not actually cluster yet: we keep them in the list clusterDigitList
 291       nextdigitC.Reset();
 292       while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) )
 293         { // scan over the list of digitsC
 294           if (digit == pMaxEnergyDigit) continue;
 295           Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
 296           AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated));
 297           //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold  )
 298           if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0  )
 299             {
 300               Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR?
 301               if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR?
 302               Bool_t shared = kFALSE; //cluster shared by 2 SuperModules?
 303               if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!!
 304                 {
 305                   clusterDigitList.AddLast(digit) ;
 306                   clusterCandidateEnergy += dEnergyCalibrated;
 307                 }
 308             }
 309         }// loop over the next digits
 310
 311       // start a cluster here only if a cluster energy is larger than clustering threshold
 312       //if (clusterCandidateEnergy > 0.1)
 313       AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold));
 314       if (clusterCandidateEnergy > fECAClusteringThreshold)
 315         {
 316           if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
 317
 318           AliEMCALRecPoint *recPoint = new  AliEMCALRecPoint("") ;
 319           fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
 320           recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
 321           fNumberOfECAClusters++ ;
 322           recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 323
 324           AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries()));
 325           for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++)
 326             {
 327
 328               digit = (AliEMCALDigit*)clusterDigitList.At(idig);
 329               Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());
 330               AliDebug(5, Form(" Adding digit %d", digit->GetId()));
 331               // note: this way the sharing info is lost!
 332               recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR?
 333               digitsC.Remove(digit);
 334             }
 335         }
 336       else
 337         {
 338           // we do not want to start clustering in the same spot!
 339           // but in this case we may NOT reuse this seed for another cluster!
 340           // need a better bookeeping?
 341           digitsC.Remove(pMaxEnergyDigit);
 342         }
 343
 344       AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries()));
 345     } // while ! done
 346
 347   //delete digitsC ; //nope we use an object
 348
 349   AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
 350 }
 351
 352 //____________________________________________________________________________
 353 void AliEMCALClusterizerNxN::MakeUnfolding()
 354 {
 355   // Unfolds clusters using the shape of an ElectroMagnetic shower
 356   // Performs unfolding of all clusters
 357
 358   if(fNumberOfECAClusters > 0){
 359     if (fGeom==0)
 360       AliFatal("Did not get geometry from EMCALLoader") ;
 361     Int_t nModulesToUnfold = fGeom->GetNCells();
 362
 363     Int_t numberofNotUnfolded = fNumberOfECAClusters ;
 364     Int_t index ;
 365     for(index = 0 ; index < numberofNotUnfolded ; index++){
 366
 367       AliEMCALRecPoint * recPoint = dynamic_cast<AliEMCALRecPoint *>( fRecPoints->At(index) ) ;
 368
 369       TVector3 gpos;
 370       Int_t absId;
 371       recPoint->GetGlobalPosition(gpos);
 372       fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId);
 373       if(absId > nModulesToUnfold)
 374         break ;
 375
 376       Int_t nMultipl = recPoint->GetMultiplicity() ;
 377       AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ;
 378       Float_t * maxAtEnergy = new Float_t[nMultipl] ;
 379       Int_t nMax = recPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ;
 380
 381       if( nMax > 1 ) {     // if cluster is very flat (no pronounced maximum) then nMax = 0
 382         //UnfoldCluster(recPoint, nMax, maxAt, maxAtEnergy) ;
 383         fRecPoints->Remove(recPoint);
 384         fRecPoints->Compress() ;
 385         index-- ;
 386         fNumberOfECAClusters-- ;
 387         numberofNotUnfolded-- ;
 388       }
 389       else{
 390         recPoint->SetNExMax(1) ; //Only one local maximum
 391       }
 392
 393       delete[] maxAt ;
 394       delete[] maxAtEnergy ;
 395     }
 396   }
 397   // End of Unfolding of clusters
 398 }
 399
 400 //____________________________________________________________________________
 401 Double_t  AliEMCALClusterizerNxN::ShowerShape(Double_t x, Double_t y)
 402 {
 403   // Shape of the shower
 404   // If you change this function, change also the gradient evaluation in ChiSquare()
 405   //
 406   Double_t r = sqrt(x*x+y*y);
 407   Double_t r133  = TMath::Power(r, 1.33) ;
 408   Double_t r669  = TMath::Power(r, 6.69) ;
 409   Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ;
 410   return shape ;
 411 }
 412
 413 //____________________________________________________________________________
 414 void  AliEMCALClusterizerNxN::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/,
 415                                             Int_t /*nMax*/,
 416                                             AliEMCALDigit ** /*maxAt*/,
 417                                             Float_t * /*maxAtEnergy*/)
 418 {
 419   //
 420   // Performs the unfolding of a cluster with nMax overlapping showers
 421   //
 422   AliWarning("Not implemented. To be.");
 423 }
 424
 425 //___________________________________________________________________
 426 void   AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val)
 427 {
 428   //
 429   // input is calibrated - the case when we run already on ESD
 430   //
 431   AliEMCALClusterizerNxN::fgkIsInputCalibrated = val;
 432 }