Allow running without tree by using
[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
483b0559 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
71
5544799a 72}
73
74//____________________________________________________________________________
40164976 75AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
ee08edde 76: AliEMCALClusterizer(geometry, calib, caloped)
0c5b726e 77{
78 // ctor, geometry and calibration are initialized elsewhere.
ee08edde 79
0c5b726e 80}
81
82
83//____________________________________________________________________________
483b0559 84 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
85{
ef305168 86 // dtor
ef305168 87}
88
89//____________________________________________________________________________
c47157cd 90void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
483b0559 91{
98e9578e 92 // Steering method to perform clusterization for the current event
93 // in AliEMCALLoader
51e4198e 94
483b0559 95 if(strstr(option,"tim"))
96 gBenchmark->Start("EMCALClusterizer");
97
98 if(strstr(option,"print"))
99 Print("") ;
51e4198e 100
1bd98442 101 //Get calibration parameters from file or digitizer default values.
102 GetCalibrationParameters() ;
51e4198e 103
40164976 104 //Get dead channel map from file or digitizer default values.
105 GetCaloCalibPedestal() ;
106
98e9578e 107 fNumberOfECAClusters = 0;
51e4198e 108
0e7c6655 109 MakeClusters() ; //only the real clusters
51e4198e 110
65bec413 111 if(fToUnfold){
112 fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
113 fClusterUnfolding->MakeUnfolding();
114 }
115
116 //Evaluate position, dispersion and other RecPoint properties for EC section
bbac7e6f 117 Int_t index ;
bbac7e6f 118 for(index = 0; index < fRecPoints->GetEntries(); index++) {
51e4198e 119 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
120 if(rp){
0d0d6b98 121 rp->EvalAll(fECAW0,fDigitsArr,fgkIsInputCalibrated) ;
51e4198e 122 //For each rec.point set the distance to the nearest bad crystal
123 rp->EvalDistanceToBadChannels(fCaloPed);
124 }
125 else AliFatal("Null rec point in list!");
bbac7e6f 126 }
51e4198e 127
bbac7e6f 128 fRecPoints->Sort() ;
51e4198e 129
bbac7e6f 130 for(index = 0; index < fRecPoints->GetEntries(); index++) {
51e4198e 131 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
132 if(rp){
133 rp->SetIndexInList(index) ;
134 rp->Print();
135 }
136 else AliFatal("Null rec point in list!");
bbac7e6f 137 }
51e4198e 138
c47157cd 139 fTreeR->Fill();
140
98e9578e 141 if(strstr(option,"deb") || strstr(option,"all"))
142 PrintRecPoints(option) ;
51e4198e 143
c47157cd 144 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
51e4198e 145
3f9ad99f 146 fRecPoints->Delete();
51e4198e 147
483b0559 148 if(strstr(option,"tim")){
149 gBenchmark->Stop("EMCALClusterizer");
98e9578e 150 printf("Exec took %f seconds for Clusterizing",
51e4198e 151 gBenchmark->GetCpuTime("EMCALClusterizer"));
98e9578e 152 }
483b0559 153}
154
155//____________________________________________________________________________
25bb3dcb 156Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
483b0559 157{
c526514b 158 // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
159 // = 1 are neighbour
160 // = 2 is in different SM; continue searching
161 // In case it is in different SM, but same phi rack, check if neigbours at eta=0
162 // neighbours are defined as digits having at least a common side
163 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
164 // which is compared to a digit (d2) not yet in a cluster
165
166 static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
167 static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
168
169 shared = kFALSE;
170
171 fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
172 fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
173 fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
174 fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
175
176 //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010
177 if(nSupMod1 != nSupMod2 ) {
178 //Check if the 2 SM are in the same PHI position (0,1), (2,3), ...
179 Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1);
180 Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2);
181
182 if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
183
184 // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2
185 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0
186 if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols;
187 else ieta2+=AliEMCALGeoParams::fgkEMCALCols;
188
189 shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
190
191 }//Different SM, same phi
192
193 Int_t rowdiff = TMath::Abs(iphi1 - iphi2);
194 Int_t coldiff = TMath::Abs(ieta1 - ieta2) ;
195
196 // neighbours with at least common side; May 11, 2007
197 if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {
198 //Diagonal?
199 //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
200
201 if (gDebug == 2)
202 printf("AliEMCALClusterizerv1::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n",
203 d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared);
204
205 return 1;
206 }//Neighbours
207 else {
208 shared = kFALSE;
209 return 2 ;
210 }//Not neighbours
483b0559 211}
212
88cb7938 213//____________________________________________________________________________
0e7c6655 214void AliEMCALClusterizerv1::MakeClusters()
483b0559 215{
216 // Steering method to construct the clusters stored in a list of Reconstructed Points
217 // A cluster is defined as a list of neighbour digits
f1487f22 218 // Mar 03, 2007 by PAI
51e4198e 219
f1487f22 220 if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
51e4198e 221
c47157cd 222 fRecPoints->Clear();
51e4198e 223
98e9578e 224 // Set up TObjArray with pointers to digits to work on
225 TObjArray *digitsC = new TObjArray();
c47157cd 226 TIter nextdigit(fDigitsArr);
98e9578e 227 AliEMCALDigit *digit;
228 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
229 digitsC->AddLast(digit);
230 }
51e4198e 231
1d46d1f6 232 double e = 0.0, ehs = 0.0;
98e9578e 233 TIter nextdigitC(digitsC);
1d46d1f6 234 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
829ba234 235 e = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId());//Time or TimeR?
236 if ( e < fMinECut) //|| digit->GetTimeR() > fTimeCut ) // time window of cell checked in calibrate
98e9578e 237 digitsC->Remove(digit);
238 else
239 ehs += e;
1d59832c 240 }
29b7e56e 241 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n",
51e4198e 242 fDigitsArr->GetEntries(),fMinECut,ehs));
243
1d46d1f6 244 nextdigitC.Reset();
51e4198e 245
1d46d1f6 246 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
c47157cd 247 TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
51e4198e 248
829ba234 249 if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()) > fECAClusteringThreshold ) ){
98e9578e 250 // start a new Tower RecPoint
c47157cd 251 if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
51e4198e 252
98e9578e 253 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
c47157cd 254 fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
255 recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
51e4198e 256 if(recPoint){
257 fNumberOfECAClusters++ ;
258
259 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
260
261 recPoint->AddDigit(*digit, Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()),kFALSE) ; //Time or TimeR?
262 TObjArray clusterDigits;
263 clusterDigits.AddLast(digit);
264 digitsC->Remove(digit) ;
265
266 AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(),
267 Calibrate(digit->GetAmplitude(),digit->GetTime(),digit->GetId()), fECAClusteringThreshold)); //Time or TimeR?
268 Float_t time = digit->GetTime();//Time or TimeR?
269 // Grow cluster by finding neighbours
270 TIter nextClusterDigit(&clusterDigits);
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
275
276 //Do not add digits with too different time
277 Bool_t shared = kFALSE;//cluster shared by 2 SuperModules?
278 if(TMath::Abs(time - digitN->GetTime()) > fTimeCut ) continue; //Time or TimeR?
279 if (AreNeighbours(digit, digitN, shared)==1) { // call (digit,digitN) in THAT order !!!!!
280 recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmplitude(), digitN->GetTime(), digitN->GetId()),shared) ;//Time or TimeR?
281 clusterDigits.AddLast(digitN) ;
282 digitsC->Remove(digitN) ;
283 } // if(ineb==1)
284 } // scan over digits
285 } // scan over digits already in cluster
286
98e9578e 287 AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
51e4198e 288 }//recpoint
289 else AliFatal("Null recpoint in array!");
98e9578e 290 } // If seed found
1d59832c 291 } // while digit
51e4198e 292
483b0559 293 delete digitsC ;
ab6a174f 294
c47157cd 295 AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));
483b0559 296}
297