EMCAL/AliEMCALClusterizerFixedWindow.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
  17
  18 // --- Root ---
  19 #include <TMath.h>
  20 #include <TMinuit.h>
  21 #include <TTree.h>
  22 #include <TBenchmark.h>
  23 #include <TBrowser.h>
  24 #include <TROOT.h>
  25 #include <TClonesArray.h>
  26 #include <TH1I.h>
  27
  28 // --- AliRoot ---
  29 #include "AliLog.h"
  30 #include "AliEMCALRecPoint.h"
  31 #include "AliEMCALDigit.h"
  32 #include "AliEMCALGeometry.h"
  33 #include "AliCaloCalibPedestal.h"
  34 #include "AliEMCALCalibData.h"
  35 #include "AliESDCaloCluster.h"
  36 #include "AliEMCALUnfolding.h"
  37
  38 #include "AliEMCALClusterizerFixedWindow.h"
  39
  40 ClassImp(AliEMCALClusterizerFixedWindow)
  41
  42 //__________________________________________________________________________________________
  43 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow() :
  44   AliEMCALClusterizer(),
  45   fNphi(4),
  46   fNeta(4),
  47   fShiftPhi(2),
  48   fShiftEta(2),
  49   fTRUshift(0),
  50   fNEtaDigitsSupMod(0),
  51   fNPhiDigitsSupMod(0),
  52   fNTRUPhi(0),
  53   fNTRUEta(0),
  54   fNEtaDigits(0),
  55   fNPhiDigits(0),
  56   fMaxShiftPhi(0),
  57   fMaxShiftEta(0),
  58   fNDigitsCluster(0),
  59   fNClusEtaNoShift(0),
  60   fNClusPhiNoShift(0),
  61   fNClusters(0),
  62   fNTotalClus(0),
  63   fClustersArray(0),
  64   fInitialized(0)
  65 {
  66   // Constructor
  67 }
  68
  69 //__________________________________________________________________________________________
  70 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry) :
  71   AliEMCALClusterizer(geometry),
  72   fNphi(4),
  73   fNeta(4),
  74   fShiftPhi(2),
  75   fShiftEta(2),
  76   fTRUshift(0),
  77   fNEtaDigitsSupMod(0),
  78   fNPhiDigitsSupMod(0),
  79   fNTRUPhi(0),
  80   fNTRUEta(0),
  81   fNEtaDigits(0),
  82   fNPhiDigits(0),
  83   fMaxShiftPhi(0),
  84   fMaxShiftEta(0),
  85   fNDigitsCluster(0),
  86   fNClusEtaNoShift(0),
  87   fNClusPhiNoShift(0),
  88   fNClusters(0),
  89   fNTotalClus(0),
  90   fClustersArray(0),
  91   fInitialized(0)
  92 {
  93   // Constructor
  94 }
  95
  96 //__________________________________________________________________________________________
  97 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped) :
  98   AliEMCALClusterizer(geometry, calib, caloped),
  99   fNphi(4),
 100   fNeta(4),
 101   fShiftPhi(2),
 102   fShiftEta(2),
 103   fTRUshift(0),
 104   fNEtaDigitsSupMod(0),
 105   fNPhiDigitsSupMod(0),
 106   fNTRUPhi(0),
 107   fNTRUEta(0),
 108   fNEtaDigits(0),
 109   fNPhiDigits(0),
 110   fMaxShiftPhi(0),
 111   fMaxShiftEta(0),
 112   fNDigitsCluster(0),
 113   fNClusEtaNoShift(0),
 114   fNClusPhiNoShift(0),
 115   fNClusters(0),
 116   fNTotalClus(0),
 117   fClustersArray(0),
 118   fInitialized(0)
 119 {
 120   // Constructor
 121 }
 122
 123 //__________________________________________________________________________________________
 124 AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
 125 {
 126   // Destructor
 127
 128   if (fClustersArray) {
 129     for (Int_t i = 0; i < fNTotalClus; i++) {
 130       if (fClustersArray[i]) {
 131         delete[] fClustersArray[i];
 132         fClustersArray[i] = 0;
 133       }
 134       delete[] fClustersArray;
 135       fClustersArray = 0;
 136     }
 137   }
 138
 139 }
 140
 141 //__________________________________________________________________________________________
 142 void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n)
 143 {
 144   // Set fNphi; if clusterizer already initialized gives a warning and does nothing
 145
 146   if (fInitialized != 0)
 147     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 148   else
 149     fNphi = n;
 150 }
 151
 152 //__________________________________________________________________________________________
 153 void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n)
 154 {
 155   // Set fNeta; if clusterizer already initialized gives a warning and does nothing
 156
 157   if (fInitialized != 0)
 158     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 159   else
 160     fNeta = n;
 161 }
 162
 163 //__________________________________________________________________________________________
 164 void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s)
 165 {
 166   // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
 167
 168   if (fInitialized != 0)
 169     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 170   else
 171     fShiftPhi = s;
 172 }
 173
 174 //__________________________________________________________________________________________
 175 void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s)
 176 {
 177   // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
 178
 179   if (fInitialized != 0)
 180     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 181   else
 182     fShiftEta = s;
 183 }
 184
 185 //__________________________________________________________________________________________
 186 void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b)
 187 {
 188   // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
 189
 190   if (fInitialized != 0)
 191     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 192   else
 193     fTRUshift = b;
 194 }
 195
 196
 197 //__________________________________________________________________________________________
 198 void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
 199 {
 200   // Steering method to perform clusterization for the current event
 201
 202   static Float_t cputime = 0;
 203   static Float_t realtime = 0;
 204
 205   if (strstr(option,"tim"))
 206     gBenchmark->Start("EMCALClusterizer");
 207
 208   if (strstr(option,"print"))
 209     Print("");
 210
 211   //Get calibration parameters from file or digitizer default values.
 212   GetCalibrationParameters();
 213
 214   //Get dead channel map from file or digitizer default values.
 215   GetCaloCalibPedestal();
 216
 217   MakeClusters();  //only the real clusters
 218
 219   if (fToUnfold) {
 220     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
 221     fClusterUnfolding->MakeUnfolding();
 222   }
 223
 224   //Evaluate position, dispersion and other RecPoint properties for EC section
 225   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
 226     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 227     if (rp) {
 228       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
 229       AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
 230       //For each rec.point set the distance to the nearest bad crystal
 231       if (fCaloPed)
 232         rp->EvalDistanceToBadChannels(fCaloPed);
 233     }
 234   }
 235
 236   fRecPoints->Sort();
 237
 238   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
 239     AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 240     if (rp) {
 241       rp->SetIndexInList(index);
 242     }
 243     else AliFatal("RecPoint NULL!!");
 244   }
 245
 246   if (fTreeR)
 247     fTreeR->Fill();
 248
 249   if (strstr(option,"deb") || strstr(option,"all"))
 250     PrintRecPoints(option);
 251
 252   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
 253
 254   if (strstr(option,"tim")) {
 255     gBenchmark->Stop("EMCALClusterizer");
 256     Printf("Exec took %f CPU time (%f real time) for clusterizing",
 257            gBenchmark->GetCpuTime("EMCALClusterizer")-cputime,gBenchmark->GetRealTime("EMCALClusterizer")-realtime);
 258     cputime = gBenchmark->GetCpuTime("EMCALClusterizer");
 259     realtime = gBenchmark->GetRealTime("EMCALClusterizer");
 260   }
 261 }
 262
 263 //__________________________________________________________________________________________
 264 void AliEMCALClusterizerFixedWindow::ExecOnce()
 265 {
 266   // Initialize clusterizer.
 267
 268   fInitialized = -1;
 269
 270   if (!fGeom) {
 271     AliError("Did not get geometry!");
 272     return;
 273   }
 274
 275   // Defining geometry and clusterization parameter
 276   fNEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
 277   fNPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
 278
 279   fNTRUPhi = 1;
 280   fNTRUEta = 1;
 281
 282   fNEtaDigits = fNEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
 283   fNPhiDigits = fNPhiDigitsSupMod * fGeom->GetNPhiSuperModule();
 284
 285   if (fTRUshift){
 286     fNTRUPhi = fGeom->GetNPhiSuperModule() * 3;
 287     fNTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
 288     fNEtaDigits /= fNTRUEta;
 289     fNPhiDigits /= fNTRUPhi;
 290   }
 291
 292   // Check if clusterizer parameter are compatible with calorimeter geometry
 293   if (fNEtaDigits < fNeta){
 294     AliError(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,fNEtaDigits));
 295     return;
 296   }
 297   if (fNPhiDigits < fNphi){
 298     AliError(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,fNPhiDigits));
 299     return;
 300   }
 301   if (fNEtaDigits % fShiftEta != 0){
 302     AliError(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,fNEtaDigits));
 303     return;
 304   }
 305   if (fNPhiDigits % fShiftPhi != 0){
 306     AliError(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,fNPhiDigits));
 307     return;
 308   }
 309   if (fNeta % fShiftEta != 0){
 310     AliError(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
 311     return;
 312   }
 313   if (fNphi % fShiftPhi != 0){
 314     AliError(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
 315     return;
 316   }
 317
 318   fMaxShiftPhi = fNphi / fShiftPhi;
 319   fMaxShiftEta = fNeta / fShiftEta;
 320
 321   fNClusEtaNoShift = fNEtaDigits / fNeta;
 322   fNClusPhiNoShift = fNPhiDigits / fNphi;
 323
 324   fNClusters = fNClusEtaNoShift * fNClusPhiNoShift * fNTRUEta * fNTRUPhi;
 325
 326   fNTotalClus = fNClusters * fMaxShiftEta * fMaxShiftPhi;
 327
 328   fNDigitsCluster = fNphi * fNeta;
 329
 330   if (fClustersArray) {
 331     for (Int_t i = 0; i < fNTotalClus; i++) {
 332       if (fClustersArray[i]) {
 333         delete[] fClustersArray[i];
 334         fClustersArray[i] = 0;
 335       }
 336       delete[] fClustersArray;
 337       fClustersArray = 0;
 338     }
 339   }
 340
 341   fClustersArray = new AliEMCALDigit**[fNTotalClus];
 342   for (Int_t i = 0; i < fNTotalClus; i++) {
 343     fClustersArray[i] = new AliEMCALDigit*[fNDigitsCluster];
 344     for (Int_t j = 0; j < fNDigitsCluster; j++) {
 345       fClustersArray[i][j] = 0;
 346     }
 347   }
 348
 349   AliDebug(1,Form("****ExecOnce*****\n"
 350                   "fNphi = %d, fNeta = %d, fShiftPhi = %d, fShiftEta = %d, fTRUshift = %d\n"
 351                   "fNEtaDigitsSupMod = %d, fNPhiDigitsSupMod = %d, fNTRUPhi = %d, fNTRUEta = %d, fNEtaDigits = %d, fNPhiDigits = %d\n"
 352                   "fMaxShiftPhi = %d, fMaxShiftEta = %d, fNDigitsCluster = %d, fNClusEtaNoShift = %d, fNClusPhiNoShift = %d\n"
 353                   "fNClusters = %d, fNTotalClus = %d\n",
 354                   fNphi,fNeta,fShiftPhi,fShiftEta,fTRUshift,
 355                   fNEtaDigitsSupMod,fNPhiDigitsSupMod,fNTRUPhi,fNTRUEta,fNEtaDigits,fNPhiDigits,
 356                   fMaxShiftPhi,fMaxShiftEta,fNDigitsCluster,fNClusEtaNoShift,fNClusPhiNoShift,
 357                   fNClusters,fNTotalClus));
 358
 359   fInitialized = 1;
 360 }
 361
 362 //__________________________________________________________________________________________
 363 void AliEMCALClusterizerFixedWindow::MakeClusters()
 364 {
 365   // Make clusters.
 366
 367   fNumberOfECAClusters = 0;
 368   fRecPoints->Delete();
 369
 370   if (fInitialized == 0)
 371     ExecOnce();
 372
 373   if (fInitialized == -1) {
 374     AliError(Form("%s: error initializing the clusterizer. No clusterization will be performed.",GetName()));
 375     return;
 376   }
 377
 378   // Set up TObjArray with pointers to digits to work on calibrated digits
 379   TObjArray *digitsC = new TObjArray();
 380   AliEMCALDigit *digit;
 381   Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
 382   TIter nextdigit(fDigitsArr);
 383   while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigit()))) { // calibrate and clean up digits
 384     dEnergyCalibrated =  digit->GetAmplitude();
 385     time              =  digit->GetTime();
 386     Calibrate(dEnergyCalibrated, time, digit->GetId());
 387     digit->SetCalibAmp(dEnergyCalibrated);
 388     digit->SetTime(time);
 389     if (dEnergyCalibrated < fMinECut || time > fTimeMax || time < fTimeMin) {
 390       continue;
 391     }
 392     else if (!fGeom->CheckAbsCellId(digit->GetId())) {
 393       continue;
 394     }
 395     else {
 396       ehs += dEnergyCalibrated;
 397       digitsC->AddLast(digit);
 398     }
 399   }
 400
 401   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
 402                   fDigitsArr->GetEntries(),fMinECut,ehs));
 403
 404   Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
 405   Int_t iphi=0, ieta=0;  // cell eta-phi indexes in SM
 406
 407   for (Int_t ishiftPhi = 0; ishiftPhi < fMaxShiftPhi; ishiftPhi++){
 408     Int_t nClusPhi = (fNPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
 409
 410     for (Int_t ishiftEta = 0; ishiftEta < fMaxShiftEta; ishiftEta++) {
 411
 412       Int_t nClusEta = (fNEtaDigits - fShiftEta * ishiftEta) / fNeta;
 413
 414       Int_t iTotalClus = fNClusters * (ishiftPhi * fMaxShiftEta + ishiftEta);
 415
 416       TIter nextdigitC(digitsC);
 417       while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
 418
 419         fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
 420         fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
 421
 422         Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + fNPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
 423
 424         Int_t iTRUphi = iphi_eff / fNPhiDigits;
 425
 426         iphi_eff -= iTRUphi * fNPhiDigits;
 427
 428         Int_t iClusPhi = iphi_eff / fNphi;
 429
 430         if (iphi_eff < 0 || iClusPhi >= nClusPhi)
 431           continue;
 432
 433         Int_t ieta_eff = ieta - fShiftEta * ishiftEta + fNEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
 434
 435         Int_t iTRUeta = ieta_eff / fNEtaDigits;
 436
 437         ieta_eff -= iTRUeta * fNEtaDigits;
 438
 439         Int_t iClusEta = ieta_eff / fNeta;
 440
 441         if (ieta_eff < 0 || iClusEta >= nClusEta)
 442           continue;
 443
 444         iphi_eff += iTRUphi * fNPhiDigits;
 445         iClusPhi = iphi_eff / fNphi;
 446
 447         ieta_eff += iTRUeta * fNEtaDigits;
 448         iClusEta = ieta_eff / fNeta;
 449
 450         Int_t iCluster = iClusPhi + iClusEta * fNClusPhiNoShift * fNTRUPhi;
 451         Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
 452
 453         if (iCluster >= fNClusters){
 454           AliError(Form("iCluster out of range! iCluster = %d, fNClusters = %d (should never happen...)", iCluster, fNClusters));
 455           return;
 456         }
 457
 458         iCluster += iTotalClus;
 459
 460         if (iCluster >= fNTotalClus){
 461           AliError(Form("iCluster out of range! iCluster = %d, fNTotalClus = %d (should never happen...)", iCluster, fNTotalClus));
 462           return;
 463         }
 464
 465         if (iDigit >= fNDigitsCluster){
 466           AliError(Form("iDigit out of range! iDigit = %d, fNDigitsCluster = %d (should never happen...)", iDigit, fNDigitsCluster));
 467           return;
 468         }
 469
 470         if (fClustersArray[iCluster][iDigit] != 0){
 471           AliError("Digit already added! (should never happen...)");
 472           return;
 473         }
 474
 475         fClustersArray[iCluster][iDigit] = digit;
 476
 477       } // loop on digit
 478
 479     } // loop on eta shift
 480
 481   } // loop on phi shift
 482
 483   AliEMCALRecPoint *recPoint = 0;
 484   Bool_t recPointOk = kFALSE;
 485   for (Int_t iCluster = 0; iCluster < fNTotalClus; iCluster++) {
 486
 487     if (!recPoint) {
 488       if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(fNumberOfECAClusters*2+1);
 489
 490       (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("");
 491       recPoint = static_cast<AliEMCALRecPoint*>(fRecPoints->At(fNumberOfECAClusters));
 492     }
 493
 494     if (recPoint) {
 495       recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 496       recPoint->SetUniqueID(iCluster);
 497       fNumberOfECAClusters++;
 498
 499       for (Int_t iDigit = 0; iDigit < fNDigitsCluster; iDigit++) {
 500         if (fClustersArray[iCluster][iDigit] == 0) continue;
 501         digit = fClustersArray[iCluster][iDigit];
 502         recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
 503         fClustersArray[iCluster][iDigit] = 0;
 504         recPointOk = kTRUE;
 505       }
 506
 507       if (recPointOk) { // unset the pointer so that a new rec point will be allocated in the next iteration
 508         recPoint = 0;
 509         recPointOk = kFALSE;
 510       }
 511     }
 512     else {
 513       AliError("Error allocating rec points!");
 514       break;
 515     }
 516   }
 517
 518   if (!recPointOk) {
 519     fRecPoints->RemoveLast();
 520     fNumberOfECAClusters--;
 521   }
 522
 523   delete digitsC;
 524   AliDebug(1, Form("MakeClusters: Number of digits %d  -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
 525   AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast()));
 526 }