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       fClustersArray(0)
  51 {
  52   // Constructor
  53 }
  54
  55 //__________________________________________________________________________________________
  56 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
  57   : AliEMCALClusterizer(geometry),
  58     fNphi(4),
  59     fNeta(4),
  60     fShiftPhi(2),
  61     fShiftEta(2),
  62     fTRUshift(0),
  63     fClustersArray(0)
  64 {
  65   // Constructor
  66 }
  67
  68 //__________________________________________________________________________________________
  69 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
  70   : AliEMCALClusterizer(geometry, calib, caloped),
  71     fNphi(4),
  72     fNeta(4),
  73     fShiftPhi(2),
  74     fShiftEta(2),
  75     fTRUshift(0),
  76     fClustersArray(0)
  77 {
  78   // Constructor
  79 }
  80
  81 //__________________________________________________________________________________________
  82 AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
  83 {
  84   // Destructor
  85
  86   delete fClustersArray;
  87 }
  88
  89 //__________________________________________________________________________________________
  90 void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n)
  91 {
  92   // Set fNphi; if clusterizer already initialized gives a warning and does nothing
  93
  94   if (fClustersArray)
  95     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
  96   else
  97     fNphi = n;
  98 }
  99
 100 //__________________________________________________________________________________________
 101 void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n)
 102 {
 103   // Set fNeta; if clusterizer already initialized gives a warning and does nothing
 104
 105   if (fClustersArray)
 106     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 107   else
 108     fNeta = n;
 109 }
 110
 111 //__________________________________________________________________________________________
 112 void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s)
 113 {
 114   // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
 115
 116   if (fClustersArray)
 117     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 118   else
 119     fShiftPhi = s;
 120 }
 121
 122 //__________________________________________________________________________________________
 123 void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s)
 124 {
 125   // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
 126
 127   if (fClustersArray)
 128     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 129   else
 130     fShiftEta = s;
 131 }
 132
 133 //__________________________________________________________________________________________
 134 void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b)
 135 {
 136   // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
 137
 138   if (fClustersArray)
 139     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
 140   else
 141     fTRUshift = b;
 142 }
 143
 144
 145 //__________________________________________________________________________________________
 146 void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
 147 {
 148   // Steering method to perform clusterization for the current event
 149
 150   if (strstr(option,"tim"))
 151     gBenchmark->Start("EMCALClusterizer");
 152
 153   if (strstr(option,"print"))
 154     Print("");
 155
 156   //Get calibration parameters from file or digitizer default values.
 157   GetCalibrationParameters();
 158
 159   //Get dead channel map from file or digitizer default values.
 160   GetCaloCalibPedestal();
 161
 162   MakeClusters();  //only the real clusters
 163
 164   if (fToUnfold) {
 165     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
 166     fClusterUnfolding->MakeUnfolding();
 167   }
 168
 169   //Evaluate position, dispersion and other RecPoint properties for EC section
 170   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
 171     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 172     if (rp) {
 173       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
 174       AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
 175       //For each rec.point set the distance to the nearest bad crystal
 176       if (fCaloPed)
 177         rp->EvalDistanceToBadChannels(fCaloPed);
 178     }
 179   }
 180
 181   fRecPoints->Sort();
 182
 183   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
 184     AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
 185     if (rp) {
 186       rp->SetIndexInList(index);
 187     }
 188     else AliFatal("RecPoint NULL!!");
 189   }
 190
 191   if (fTreeR)
 192     fTreeR->Fill();
 193
 194   if (strstr(option,"deb") || strstr(option,"all"))
 195     PrintRecPoints(option);
 196
 197   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
 198
 199   if (strstr(option,"tim")) {
 200     gBenchmark->Stop("EMCALClusterizer");
 201     printf("Exec took %f seconds for Clusterizing",
 202            gBenchmark->GetCpuTime("EMCALClusterizer"));
 203   }
 204 }
 205
 206 //__________________________________________________________________________________________
 207 void AliEMCALClusterizerFixedWindow::MakeClusters()
 208 {
 209   // Make clusters
 210
 211   if (fGeom == 0)
 212     AliFatal("Did not get geometry from EMCALLoader");
 213
 214   fNumberOfECAClusters = 0;
 215   fRecPoints->Delete();
 216
 217   Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphi=0, ieta=0;
 218
 219   // Defining geometry and clusterization parameter
 220   Int_t nEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
 221   Int_t nPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
 222
 223   Int_t nTRUPhi = 1;
 224   Int_t nTRUEta = 1;
 225
 226   Int_t nEtaDigits = nEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
 227   Int_t nPhiDigits = nPhiDigitsSupMod * fGeom->GetNPhiSuperModule();
 228
 229   if (fTRUshift){
 230     nTRUPhi = fGeom->GetNPhiSuperModule() * 3;
 231     nTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
 232     nEtaDigits /= nTRUEta;
 233     nPhiDigits /= nTRUPhi;
 234   }
 235
 236   // Check if clusterizer parameter are compatible with calorimeter geometry
 237   if (nEtaDigits < fNeta){
 238     AliFatal(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,nEtaDigits));
 239     return;
 240   }
 241   if (nPhiDigits < fNphi){
 242     AliFatal(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,nPhiDigits));
 243     return;
 244   }
 245   if (nEtaDigits % fShiftEta != 0){
 246     AliFatal(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,nEtaDigits));
 247     return;
 248   }
 249   if (nPhiDigits % fShiftPhi != 0){
 250     AliFatal(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,nPhiDigits));
 251     return;
 252   }
 253   if (fNeta % fShiftEta != 0){
 254     AliFatal(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
 255     return;
 256   }
 257   if (fNphi % fShiftPhi != 0){
 258     AliFatal(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
 259     return;
 260   }
 261
 262   Int_t maxiShiftPhi = fNphi / fShiftPhi;
 263   Int_t maxiShiftEta = fNeta / fShiftEta;
 264
 265   Int_t nDigitsCluster = fNphi * fNeta;
 266
 267   Int_t nClusEtaNoShift = nEtaDigits / fNeta;
 268   Int_t nClusPhiNoShift = nPhiDigits / fNphi;
 269
 270   Int_t nClusters =  nClusEtaNoShift * nClusPhiNoShift * nTRUEta * nTRUPhi;
 271
 272   Int_t nTotalClus = nClusters * maxiShiftEta * maxiShiftPhi;
 273
 274   if (!fClustersArray) {
 275     fClustersArray = new AliEMCALDigit**[nTotalClus];
 276     for (Int_t i = 0; i < nTotalClus; i++)
 277     {
 278       fClustersArray[i] = NULL;
 279     }
 280   }
 281
 282   // Set up TObjArray with pointers to digits to work on calibrated digits
 283   TObjArray *digitsC = new TObjArray();
 284   AliEMCALDigit *digit;
 285   Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
 286   TIter nextdigit(fDigitsArr);
 287   while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigit()))) { // calibrate and clean up digits
 288     dEnergyCalibrated =  digit->GetAmplitude();
 289     time              =  digit->GetTime();
 290     Calibrate(dEnergyCalibrated, time, digit->GetId());
 291     digit->SetCalibAmp(dEnergyCalibrated);
 292     digit->SetTime(time);
 293     if (dEnergyCalibrated < fMinECut || time > fTimeMax || time < fTimeMin) {
 294       continue;
 295     }
 296     else if (!fGeom->CheckAbsCellId(digit->GetId())) {
 297       continue;
 298     }
 299     else {
 300       ehs += dEnergyCalibrated;
 301       digitsC->AddLast(digit);
 302     }
 303   }
 304
 305   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
 306                   fDigitsArr->GetEntries(),fMinECut,ehs));
 307
 308   for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++){
 309     Int_t nClusPhi = (nPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
 310
 311     for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++) {
 312
 313       Int_t nClusEta = (nEtaDigits - fShiftEta * ishiftEta) / fNeta;
 314
 315       Int_t iTotalClus = nClusters * (ishiftPhi * maxiShiftEta + ishiftEta);
 316
 317       TIter nextdigitC(digitsC);
 318       while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
 319
 320         fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
 321         fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
 322
 323         Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
 324
 325         Int_t iTRUphi = iphi_eff / nPhiDigits;
 326
 327         iphi_eff -= iTRUphi * nPhiDigits;
 328
 329         Int_t iClusPhi = iphi_eff / fNphi;
 330
 331         if (iphi_eff < 0 || iClusPhi >= nClusPhi)
 332           continue;
 333
 334         Int_t ieta_eff = ieta - fShiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
 335
 336         Int_t iTRUeta = ieta_eff / nEtaDigits;
 337
 338         ieta_eff -= iTRUeta * nEtaDigits;
 339
 340         Int_t iClusEta = ieta_eff / fNeta;
 341
 342         if (ieta_eff < 0 || iClusEta >= nClusEta)
 343           continue;
 344
 345         iphi_eff += iTRUphi * nPhiDigits;
 346         iClusPhi = iphi_eff / fNphi;
 347
 348         ieta_eff += iTRUeta * nEtaDigits;
 349         iClusEta = ieta_eff / fNeta;
 350
 351         Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi;
 352         Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
 353
 354         if (iCluster >= nClusters){
 355           AliWarning(Form("iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
 356           return;
 357         }
 358
 359         iCluster += iTotalClus;
 360
 361         if (fClustersArray[iCluster] == NULL){
 362           fNumberOfECAClusters++;
 363           fClustersArray[iCluster] = new AliEMCALDigit*[nDigitsCluster];
 364           for (Int_t i = 0; i < nDigitsCluster; i++){
 365             fClustersArray[iCluster][i] = NULL;
 366           }
 367         }
 368
 369         if (fClustersArray[iCluster][iDigit] != NULL){
 370           AliWarning("Digit already added!");
 371           return;
 372         }
 373
 374         fClustersArray[iCluster][iDigit] = digit;
 375
 376       } // loop on digit
 377
 378     } // loop on eta shift
 379
 380   } // loop on phi shift
 381
 382   if(fNumberOfECAClusters >= fRecPoints->GetSize())
 383     fRecPoints->Expand(fNumberOfECAClusters+1);
 384
 385   Int_t iRecPoint = 0;
 386   for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++) {
 387
 388     if (fClustersArray[iCluster] == NULL)
 389       continue;
 390
 391     (*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
 392     AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint*> (fRecPoints->At(iRecPoint));
 393
 394     if (recPoint) {
 395       iRecPoint++;
 396       recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
 397       recPoint->SetUniqueID(iCluster);
 398
 399       for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++){
 400         if (fClustersArray[iCluster][iDigit] == NULL) continue;
 401         digit = fClustersArray[iCluster][iDigit];
 402         recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
 403         fClustersArray[iCluster][iDigit] = NULL;
 404       }
 405     }
 406
 407     delete[] fClustersArray[iCluster];
 408     fClustersArray[iCluster] = NULL;
 409   }
 410
 411   delete digitsC;
 412   AliDebug(1, Form("MakeClusters: Number of digits %d  -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
 413   AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast()));
 414 }