Adding include path to allow compilation of CleanGeom task
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerv1.cxx
CommitLineData
483b0559 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 **************************************************************************/
173558f2 15
483b0559 16/* $Id$ */
803d1ab0 17
3a8be91c 18//-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
ee08edde 19//-- Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
483b0559 20//////////////////////////////////////////////////////////////////////////////
21// Clusterization class. Performs clusterization (collects neighbouring active cells) and
22// unfolds the clusters having several local maxima.
23// Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
24// EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
25// parameters including input digits branch title, thresholds etc.)
ee08edde 26//
483b0559 27
28// --- ROOT system ---
29
e52475ed 30#include <TFile.h>
e52475ed 31#include <TMath.h>
32#include <TMinuit.h>
33#include <TTree.h>
e52475ed 34#include <TBenchmark.h>
35#include <TBrowser.h>
a1e17193 36#include <TROOT.h>
ee08edde 37#include <TList.h>
38#include <TClonesArray.h>
1d59832c 39
483b0559 40// --- Standard library ---
ee08edde 41#include <cassert>
173558f2 42
483b0559 43// --- AliRoot header files ---
ee08edde 44#include "AliLog.h"
483b0559 45#include "AliEMCALClusterizerv1.h"
70a93198 46#include "AliEMCALRecPoint.h"
483b0559 47#include "AliEMCALDigit.h"
05a92d59 48#include "AliEMCALGeometry.h"
40164976 49#include "AliCaloCalibPedestal.h"
a435f763 50#include "AliEMCALCalibData.h"
ee08edde 51#include "AliESDCaloCluster.h"
65bec413 52#include "AliEMCALUnfolding.h"
483b0559 53
54ClassImp(AliEMCALClusterizerv1)
1963b290 55
483b0559 56//____________________________________________________________________________
ee08edde 57AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer()
483b0559 58{
59 // ctor with the indication of the file where header Tree and digits Tree are stored
60
62c48b86 61 Init();
483b0559 62}
05a92d59 63
483b0559 64//____________________________________________________________________________
5544799a 65AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
ee08edde 66 : AliEMCALClusterizer(geometry)
5544799a 67{
68 // ctor with the indication of the file where header Tree and digits Tree are stored
69 // use this contructor to avoid usage of Init() which uses runloader
70 // change needed by HLT - MP
5544799a 71}
72
73//____________________________________________________________________________
40164976 74AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
ee08edde 75: AliEMCALClusterizer(geometry, calib, caloped)
0c5b726e 76{
62c48b86 77 // ctor, geometry and calibration are initialized elsewhere.
0c5b726e 78}
79
0c5b726e 80//____________________________________________________________________________
483b0559 81 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
82{
ef305168 83 // dtor
ef305168 84}
85
86//____________________________________________________________________________
c47157cd 87void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
483b0559 88{
98e9578e 89 // Steering method to perform clusterization for the current event
90 // in AliEMCALLoader
51e4198e 91
483b0559 92 if(strstr(option,"tim"))
93 gBenchmark->Start("EMCALClusterizer");
94
95 if(strstr(option,"print"))
62c48b86 96 Print("");
51e4198e 97
1bd98442 98 //Get calibration parameters from file or digitizer default values.
62c48b86 99 GetCalibrationParameters();
51e4198e 100
40164976 101 //Get dead channel map from file or digitizer default values.
62c48b86 102 GetCaloCalibPedestal();
40164976 103
98e9578e 104 fNumberOfECAClusters = 0;
51e4198e 105
62c48b86 106 MakeClusters(); //only the real clusters
51e4198e 107
65bec413 108 if(fToUnfold){
109 fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
110 fClusterUnfolding->MakeUnfolding();
111 }
112
113 //Evaluate position, dispersion and other RecPoint properties for EC section
62c48b86 114 Int_t index;
bbac7e6f 115 for(index = 0; index < fRecPoints->GetEntries(); index++) {
51e4198e 116 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
117 if(rp){
294553d1 118 rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
51e4198e 119 //For each rec.point set the distance to the nearest bad crystal
62c48b86 120 if (fCaloPed)
121 rp->EvalDistanceToBadChannels(fCaloPed);
51e4198e 122 }
123 else AliFatal("Null rec point in list!");
bbac7e6f 124 }
51e4198e 125
62c48b86 126 fRecPoints->Sort();
51e4198e 127
bbac7e6f 128 for(index = 0; index < fRecPoints->GetEntries(); index++) {
51e4198e 129 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
130 if(rp){
62c48b86 131 rp->SetIndexInList(index);
51e4198e 132 rp->Print();
133 }
134 else AliFatal("Null rec point in list!");
bbac7e6f 135 }
51e4198e 136
62c48b86 137 if (fTreeR)
138 fTreeR->Fill();
c47157cd 139
98e9578e 140 if(strstr(option,"deb") || strstr(option,"all"))
62c48b86 141 PrintRecPoints(option);
51e4198e 142
c47157cd 143 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
51e4198e 144
483b0559 145 if(strstr(option,"tim")){
146 gBenchmark->Stop("EMCALClusterizer");
98e9578e 147 printf("Exec took %f seconds for Clusterizing",
51e4198e 148 gBenchmark->GetCpuTime("EMCALClusterizer"));
98e9578e 149 }
483b0559 150}
151
152//____________________________________________________________________________
25bb3dcb 153Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
483b0559 154{
62c48b86 155 // Gives the neighbourness of two digits = 0 are not neighbour; continue searching
c526514b 156 // = 1 are neighbour
157 // = 2 is in different SM; continue searching
158 // In case it is in different SM, but same phi rack, check if neigbours at eta=0
159 // neighbours are defined as digits having at least a common side
160 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
161 // which is compared to a digit (d2) not yet in a cluster
162
62c48b86 163 Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
164 Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
c526514b 165
166 shared = kFALSE;
167
168 fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
169 fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
170 fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
171 fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
172
173 //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
62c48b86 174 if (nSupMod1 != nSupMod2 ) {
c526514b 175 //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
176 Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
177 Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
178
179 if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
180
181 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
182 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
183 if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
184 else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
185
186 shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
62c48b86 187 } //Different SM, same phi
c526514b 188
189 Int_t rowdiff = TMath::Abs(iphi1 - iphi2);
62c48b86 190 Int_t coldiff = TMath::Abs(ieta1 - ieta2);
c526514b 191
192 // neighbours with at least common side; May 11, 2007
193 if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {
194 //Diagonal?
195 //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
196
197 if (gDebug == 2)
198 printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
199 d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);
c526514b 200 return 1;
62c48b86 201 } //Neighbours
c526514b 202 else {
203 shared = kFALSE;
62c48b86 204 return 2;
205 } //Not neighbours
483b0559 206}
207
88cb7938 208//____________________________________________________________________________
0e7c6655 209void AliEMCALClusterizerv1::MakeClusters()
483b0559 210{
211 // Steering method to construct the clusters stored in a list of Reconstructed Points
212 // A cluster is defined as a list of neighbour digits
f1487f22 213 // Mar 03, 2007 by PAI
51e4198e 214
f1487f22 215 if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
51e4198e 216
62c48b86 217 fRecPoints->Delete();
51e4198e 218
783153ff 219 // Set up TObjArray with pointers to digits to work on calibrated digits
98e9578e 220 TObjArray *digitsC = new TObjArray();
98e9578e 221 AliEMCALDigit *digit;
783153ff 222 Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
223 TIter nextdigit(fDigitsArr);
224 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // calibrate and clean up digits
225 dEnergyCalibrated = digit->GetAmplitude();
226 time = digit->GetTime();
227 Calibrate(dEnergyCalibrated,time,digit->GetId());
228 digit->SetCalibAmp(dEnergyCalibrated);
229 digit->SetTime(time);
959d5ca4 230
231 if ( dEnergyCalibrated < fMinECut || time > fTimeMax || time < fTimeMin ){
783153ff 232 continue;
233 }
234 else if (!fGeom->CheckAbsCellId(digit->GetId()))
235 continue;
236 else{
237 ehs += dEnergyCalibrated;
238 digitsC->AddLast(digit);
239 }
1d59832c 240 }
783153ff 241
29b7e56e 242 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n",
51e4198e 243 fDigitsArr->GetEntries(),fMinECut,ehs));
244
783153ff 245 TIter nextdigitC(digitsC);
1d46d1f6 246 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
c47157cd 247 TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
783153ff 248 dEnergyCalibrated = digit->GetCalibAmp();
249 time = digit->GetTime();
250 if(fGeom->CheckAbsCellId(digit->GetId()) && ( dEnergyCalibrated > fECAClusteringThreshold ) ){
98e9578e 251 // start a new Tower RecPoint
62c48b86 252 if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1);
51e4198e 253
62c48b86 254 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
255 fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
256 recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters));
257 if (recPoint) {
258 fNumberOfECAClusters++;
51e4198e 259
260 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
783153ff 261 recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
51e4198e 262 TObjArray clusterDigits;
263 clusterDigits.AddLast(digit);
62c48b86 264 digitsC->Remove(digit);
51e4198e 265
783153ff 266 AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(), dEnergyCalibrated, fECAClusteringThreshold)); //Time or TimeR?
267
51e4198e 268 // Grow cluster by finding neighbours
269 TIter nextClusterDigit(&clusterDigits);
783153ff 270
51e4198e 271 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
272 TIter nextdigitN(digitsC);
273 AliEMCALDigit *digitN = 0; // digi neighbor
274 while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
51e4198e 275 //Do not add digits with too different time
276 Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
277 if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
278 if (AreNeighbours(digit, digitN, shared)==1) { // call (digit,digitN) in THAT order !!!!!
783153ff 279 recPoint->AddDigit(*digitN, digitN->GetCalibAmp(), shared);
62c48b86 280 clusterDigits.AddLast(digitN);
281 digitsC->Remove(digitN);
51e4198e 282 } // if(ineb==1)
283 } // scan over digits
284 } // scan over digits already in cluster
285
98e9578e 286 AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
51e4198e 287 }//recpoint
288 else AliFatal("Null recpoint in array!");
98e9578e 289 } // If seed found
1d59832c 290 } // while digit
51e4198e 291
62c48b86 292 delete digitsC;
ab6a174f 293
c47157cd 294 AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
483b0559 295}