--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+// This class derives from AliEMCALClustrerizer but also keeps the API of AliEMCALClusterizerv1
+// Algorithm:
+// 1. peek the most energetic cell
+// 2. assign it as a center of the cluster and add cells surrounding it: 3x3, 5x5...
+// 3. remove the cells contributing to the cluster
+// 4. start from 1 for the remaining clusters
+// 5. cluster splitting (not implemented yet) - use the shape analysis to resolve the energy sharing
+// - for high energy clusters check the surrounding of the 3x3 clusters for extra energy
+// (merge 3x3 clusters and resolve the internal energy sharing - case for 2 clusters merged)
+// Use Case:
+// root [0] AliEMCALClusterizerFixedWindow * cl = new AliEMCALClusterizerFixedWindow("galice.root")
+// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
+// //reads gAlice from header file "..."
+// root [1] cl->ExecuteTask()
+// //finds RecPoints in all events stored in galice.root
+// root [2] cl->SetDigitsBranch("digits2")
+// //sets another title for Digitis (input) branch
+// root [3] cl->SetRecPointsBranch("recp2")
+// //sets another title four output branches
+// root [4] cl->SetTowerLocalMaxCut(0.03)
+// //set clusterization parameters
+// root [5] cl->ExecuteTask("deb all time")
+// //once more finds RecPoints options are
+// // deb - print number of found rec points
+// // deb all - print number of found RecPoints and some their characteristics
+// // time - print benchmarking results
+
+#include "AliEMCALClusterizerFixedWindow.h"
+
+// --- ROOT system ---
+#include <TMath.h>
+#include <TMinuit.h>
+#include <TTree.h>
+#include <TBenchmark.h>
+#include <TBrowser.h>
+#include <TROOT.h>
+#include <TClonesArray.h>
+#include <TH1I.h>
+
+// --- Standard library ---
+#include <cassert>
+//#include <iostream>
+//#include <fstream>
+
+// --- AliRoot header files ---
+#include "AliLog.h"
+#include "AliEMCALRecPoint.h"
+#include "AliEMCALDigit.h"
+#include "AliEMCALGeometry.h"
+#include "AliCaloCalibPedestal.h"
+#include "AliEMCALCalibData.h"
+#include "AliESDCaloCluster.h"
+#include "AliEMCALUnfolding.h"
+#include "AliEMCALFixedWindowClusterInfo.h"
+
+ClassImp(AliEMCALClusterizerFixedWindow)
+
+//____________________________________________________________________________
+AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow()
+: AliEMCALClusterizer(),
+nPhi(4),
+nEta(4),
+shiftPhi(2),
+shiftEta(2),
+fTRUshift(0),
+clusters_array(0),
+fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
+{
+ // ctor with the indication of the file where header Tree and digits Tree are stored
+
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
+: AliEMCALClusterizer(geometry),
+nPhi(4),
+nEta(4),
+shiftPhi(2),
+shiftEta(2),
+fTRUshift(0),
+clusters_array(0),
+fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
+{
+ // ctor with the indication of the file where header Tree and digits Tree are stored
+ // use this contructor to avoid usage of Init() which uses runloader
+ // change needed by HLT - MP
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
+: AliEMCALClusterizer(geometry, calib, caloped),
+nPhi(4),
+nEta(4),
+shiftPhi(2),
+shiftEta(2),
+fTRUshift(0),
+clusters_array(0),
+fClustersInfo(new AliEMCALFixedWindowClusterInfo("clustersInfo"))
+{
+ // ctor, geometry and calibration are initialized elsewhere.
+}
+
+//____________________________________________________________________________
+AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
+{
+ // dtor
+ delete fClustersInfo;
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
+{
+ // Steering method to perform clusterization for the current event
+ // in AliEMCALLoader
+
+ if (strstr(option,"tim"))
+ gBenchmark->Start("EMCALClusterizer");
+
+ if (strstr(option,"print"))
+ Print("");
+
+ //Get calibration parameters from file or digitizer default values.
+ GetCalibrationParameters();
+
+ //Get dead channel map from file or digitizer default values.
+ GetCaloCalibPedestal();
+
+ MakeClusters(); //only the real clusters
+
+ if (fToUnfold) {
+ fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
+ fClusterUnfolding->MakeUnfolding();
+ }
+
+ //Evaluate position, dispersion and other RecPoint properties for EC section
+ for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
+ AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
+ if (rp) {
+ rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
+ AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
+ //For each rec.point set the distance to the nearest bad crystal
+ if (fCaloPed)
+ rp->EvalDistanceToBadChannels(fCaloPed);
+ }
+ }
+
+ //fRecPoints->Sort();
+
+ for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
+ AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
+ if (rp) {
+ rp->SetIndexInList(index);
+ }
+ else AliFatal("RecPoint NULL!!");
+ }
+
+ if (fTreeR)
+ fTreeR->Fill();
+
+ if (strstr(option,"deb") || strstr(option,"all"))
+ PrintRecPoints(option);
+
+ AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
+
+ if (strstr(option,"tim")) {
+ gBenchmark->Stop("EMCALClusterizer");
+ printf("Exec took %f seconds for Clusterizing",
+ gBenchmark->GetCpuTime("EMCALClusterizer"));
+ }
+}
+
+//____________________________________________________________________________
+void AliEMCALClusterizerFixedWindow::MakeClusters()
+{
+ // Make clusters
+
+ if (fGeom == 0)
+ AliFatal("Did not get geometry from EMCALLoader");
+
+ fNumberOfECAClusters = 0;
+ fRecPoints->Delete();
+
+ if (fClustersInfo->GetLastElementId() > 0)
+ fClustersInfo->Clear();
+
+ Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphi=0, ieta=0;
+
+ // Defining geometry and clusterization parameter
+ Int_t nEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
+ Int_t nPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
+
+ Int_t nTRUPhi = 1;
+ Int_t nTRUEta = 1;
+
+ Int_t nEtaDigits = nEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
+ Int_t nPhiDigits = nPhiDigitsSupMod * fGeom->GetNPhiSuperModule();
+
+ if (fTRUshift)
+ {
+ nTRUPhi = fGeom->GetNPhiSuperModule() * 3;
+ nTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
+ nEtaDigits /= nTRUEta;
+ nPhiDigits /= nTRUPhi;
+ }
+
+ // Check if clusterizer parameter are compatible with calorimeter geometry
+ if (nEtaDigits < nEta)
+ {
+ AliFatal(Form("Error: nEta = %d is greater than nEtaDigits = %d.",nEta,nEtaDigits));
+ return;
+ }
+ if (nPhiDigits < nPhi)
+ {
+ AliFatal(Form("Error: nPhi = %d is greater than nPhiDigits = %d.",nPhi,nPhiDigits));
+ return;
+ }
+ if (nEtaDigits % shiftEta != 0)
+ {
+ AliFatal(Form("Error: shiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",shiftEta,nEtaDigits));
+ return;
+ }
+ if (nPhiDigits % shiftPhi != 0)
+ {
+ AliFatal(Form("Error: shiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",shiftPhi,nPhiDigits));
+ return;
+ }
+ if (nEta % shiftEta != 0)
+ {
+ AliFatal(Form("Error: shiftEta = %d is not divisor of nEta = %d.",shiftEta,nEta));
+ return;
+ }
+ if (nPhi % shiftPhi != 0)
+ {
+ AliFatal(Form("Error: shiftPhi = %d is not divisor of nPhi = %d).",shiftPhi,nPhi));
+ return;
+ }
+
+ Int_t maxiShiftPhi = nPhi / shiftPhi;
+ Int_t maxiShiftEta = nEta / shiftEta;
+
+ Int_t nDigitsCluster = nPhi * nEta;
+
+ Int_t nClusEtaNoShift = nEtaDigits / nEta;
+ Int_t nClusPhiNoShift = nPhiDigits / nPhi;
+
+ Int_t nClusters = nClusEtaNoShift * nClusPhiNoShift * nTRUEta * nTRUPhi;
+
+ Int_t nTotalClus = nClusters * maxiShiftEta * maxiShiftPhi;
+
+ if (!clusters_array)
+ {
+ clusters_array = new AliEMCALDigit**[nTotalClus];
+ for (Int_t i = 0; i < nTotalClus; i++)
+ {
+ clusters_array[i] = NULL;
+ }
+ }
+
+ AliEMCALDigit *digit = 0;
+
+ for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++)
+ {
+ Int_t nClusPhi = (nPhiDigits - shiftPhi * ishiftPhi) / nPhi;
+
+ for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++)
+ {
+
+ Int_t nClusEta = (nEtaDigits - shiftEta * ishiftEta) / nEta;
+
+ Int_t iTotalClus = nClusters * (ishiftPhi * maxiShiftEta + ishiftEta);
+
+ TIter nextdigit(fDigitsArr);
+
+ nextdigit.Reset();
+
+ while (digit = static_cast<AliEMCALDigit*>(nextdigit()))
+ {
+ fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
+ fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
+
+ Int_t iphi_eff = iphi - shiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
+
+ Int_t iTRUphi = iphi_eff / nPhiDigits;
+
+ iphi_eff -= iTRUphi * nPhiDigits;
+
+ Int_t iClusPhi = iphi_eff / nPhi;
+
+ if (iphi_eff < 0 || iClusPhi >= nClusPhi)
+ continue;
+
+ Int_t ieta_eff = ieta - shiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
+
+ Int_t iTRUeta = ieta_eff / nEtaDigits;
+
+ ieta_eff -= iTRUeta * nEtaDigits;
+
+ Int_t iClusEta = ieta_eff / nEta;
+
+ if (ieta_eff < 0 || iClusEta >= nClusEta)
+ continue;
+
+ iphi_eff += iTRUphi * nPhiDigits;
+ iClusPhi = iphi_eff / nPhi;
+
+ ieta_eff += iTRUeta * nEtaDigits;
+ iClusEta = ieta_eff / nEta;
+
+ Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi;
+ Int_t iDigit = iphi_eff % nPhi + (ieta_eff % nEta) * nPhi;
+
+
+ if (iCluster >= nClusters)
+ {
+ AliFatal(Form("ERROR: iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
+ return;
+ }
+
+ iCluster += iTotalClus;
+
+ if (clusters_array[iCluster] == NULL)
+ {
+ fNumberOfECAClusters++;
+ clusters_array[iCluster] = new AliEMCALDigit*[nDigitsCluster];
+ for (Int_t i = 0; i < nDigitsCluster; i++)
+ {
+ clusters_array[iCluster][i] = NULL;
+ }
+
+ fClustersInfo->Add(iCluster, -1, iClusEta, iClusPhi);
+ }
+
+ if (clusters_array[iCluster][iDigit] != NULL)
+ {
+ AliFatal("ERROR: digit already added!");
+ return;
+ }
+
+ clusters_array[iCluster][iDigit] = digit;
+
+ } // loop on digit
+
+ } // loop on eta shift
+
+ } // loop on phi shift
+
+ Int_t iRecPoint = 0;
+ for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++)
+ {
+
+ if (clusters_array[iCluster] == NULL) continue;
+
+ (*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
+ AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint*> (fRecPoints->At(iRecPoint));
+
+ if (recPoint)
+ {
+ if (fClustersInfo->ContainsIndex(iRecPoint))
+ AliFatal(Form("ERROR: index present already, %d", iRecPoint));
+
+ fClustersInfo->SetIndexFromId(iCluster, iRecPoint);
+
+ iRecPoint++;
+ recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
+ // note: this way the sharing info is lost!
+ for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++)
+ {
+ if (clusters_array[iCluster][iDigit] == NULL) continue;
+ digit = clusters_array[iCluster][iDigit];
+ Float_t dEnergyCalibrated = digit->GetAmplitude();
+ Float_t time = digit->GetTime();
+ Calibrate(dEnergyCalibrated,time,digit->GetId());
+ digit->SetCalibAmp(dEnergyCalibrated);
+ recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE); //Time or TimeR?
+ clusters_array[iCluster][iDigit] = NULL;
+ }
+ }
+
+ delete[] clusters_array[iCluster];
+ clusters_array[iCluster] = NULL;
+ }
+
+ AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f)\n",
+ fDigitsArr->GetEntries(),fMinECut));
+
+ AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff