bug #84425: printf statements in AliTPCcalibDButil.cxx
[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
98e9578e 219 // Set up TObjArray with pointers to digits to work on
220 TObjArray *digitsC = new TObjArray();
c47157cd 221 TIter nextdigit(fDigitsArr);
98e9578e 222 AliEMCALDigit *digit;
223 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
224 digitsC->AddLast(digit);
225 }
51e4198e 226
1d46d1f6 227 double e = 0.0, ehs = 0.0;
98e9578e 228 TIter nextdigitC(digitsC);
1d46d1f6 229 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
829ba234 230 e = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());//Time or TimeR?
231 if ( e < fMinECut) //|| digit->GetTimeR() > fTimeCut ) // time window of cell checked in calibrate
98e9578e 232 digitsC->Remove(digit);
233 else
234 ehs += e;
1d59832c 235 }
29b7e56e 236 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n",
51e4198e 237 fDigitsArr->GetEntries(),fMinECut,ehs));
238
1d46d1f6 239 nextdigitC.Reset();
51e4198e 240
1d46d1f6 241 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
c47157cd 242 TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
51e4198e 243
829ba234 244 if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()) > fECAClusteringThreshold ) ){
98e9578e 245 // start a new Tower RecPoint
62c48b86 246 if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1);
51e4198e 247
62c48b86 248 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("");
249 fRecPoints->AddAt(recPoint, fNumberOfECAClusters);
250 recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters));
251 if (recPoint) {
252 fNumberOfECAClusters++;
51e4198e 253
254 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
255
62c48b86 256 recPoint->AddDigit(*digit, Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()),kFALSE); //Time or TimeR?
51e4198e 257 TObjArray clusterDigits;
258 clusterDigits.AddLast(digit);
62c48b86 259 digitsC->Remove(digit);
51e4198e 260
261 AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(),
262 Calibrate(digit->GetAmplitude(),digit->GetTime(),digit->GetId()), fECAClusteringThreshold)); //Time or TimeR?
263 Float_t time = digit->GetTime();//Time or TimeR?
264 // Grow cluster by finding neighbours
265 TIter nextClusterDigit(&clusterDigits);
266 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
267 TIter nextdigitN(digitsC);
268 AliEMCALDigit *digitN = 0; // digi neighbor
269 while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
270
271 //Do not add digits with too different time
272 Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
273 if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
274 if (AreNeighbours(digit, digitN, shared)==1) { // call (digit,digitN) in THAT order !!!!!
62c48b86 275 recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmplitude(), digitN->GetTime(), digitN->GetId()),shared); //Time or TimeR?
276 clusterDigits.AddLast(digitN);
277 digitsC->Remove(digitN);
51e4198e 278 } // if(ineb==1)
279 } // scan over digits
280 } // scan over digits already in cluster
281
98e9578e 282 AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
51e4198e 283 }//recpoint
284 else AliFatal("Null recpoint in array!");
98e9578e 285 } // If seed found
1d59832c 286 } // while digit
51e4198e 287
62c48b86 288 delete digitsC;
ab6a174f 289
c47157cd 290 AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
483b0559 291}