1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
19 //-- Gustavo Conesa (LPSC-Grenoble), move common clusterizer functionalities to mother class
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.)
28 // --- ROOT system ---
34 #include <TBenchmark.h>
38 #include <TClonesArray.h>
40 // --- Standard library ---
43 // --- AliRoot header files ---
45 #include "AliEMCALClusterizerv1.h"
46 #include "AliEMCALRecPoint.h"
47 #include "AliEMCALDigit.h"
48 #include "AliEMCALGeometry.h"
49 #include "AliCaloCalibPedestal.h"
50 #include "AliEMCALCalibData.h"
51 #include "AliESDCaloCluster.h"
52 #include "AliEMCALUnfolding.h"
54 ClassImp(AliEMCALClusterizerv1)
56 //____________________________________________________________________________
57 AliEMCALClusterizerv1::AliEMCALClusterizerv1(): AliEMCALClusterizer()
59 // ctor with the indication of the file where header Tree and digits Tree are stored
64 //____________________________________________________________________________
65 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry)
66 : AliEMCALClusterizer(geometry)
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
74 //____________________________________________________________________________
75 AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
76 : AliEMCALClusterizer(geometry, calib, caloped)
78 // ctor, geometry and calibration are initialized elsewhere.
83 //____________________________________________________________________________
84 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
89 //____________________________________________________________________________
90 void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option)
92 // Steering method to perform clusterization for the current event
95 if(strstr(option,"tim"))
96 gBenchmark->Start("EMCALClusterizer");
98 if(strstr(option,"print"))
101 //Get calibration parameters from file or digitizer default values.
102 GetCalibrationParameters() ;
104 //Get dead channel map from file or digitizer default values.
105 GetCaloCalibPedestal() ;
107 fNumberOfECAClusters = 0;
109 MakeClusters() ; //only the real clusters
112 fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
113 fClusterUnfolding->MakeUnfolding();
116 //Evaluate position, dispersion and other RecPoint properties for EC section
118 for(index = 0; index < fRecPoints->GetEntries(); index++) {
119 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
121 rp->EvalAll(fECAW0,fDigitsArr) ;
122 //For each rec.point set the distance to the nearest bad crystal
123 rp->EvalDistanceToBadChannels(fCaloPed);
125 else AliFatal("Null rec point in list!");
130 for(index = 0; index < fRecPoints->GetEntries(); index++) {
131 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
133 rp->SetIndexInList(index) ;
136 else AliFatal("Null rec point in list!");
141 if(strstr(option,"deb") || strstr(option,"all"))
142 PrintRecPoints(option) ;
144 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
146 fRecPoints->Delete();
148 if(strstr(option,"tim")){
149 gBenchmark->Stop("EMCALClusterizer");
150 printf("Exec took %f seconds for Clusterizing",
151 gBenchmark->GetCpuTime("EMCALClusterizer"));
155 //____________________________________________________________________________
156 Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const
158 // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
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
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;
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);
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);
182 if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours
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;
189 shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment.
191 }//Different SM, same phi
193 Int_t rowdiff = TMath::Abs(iphi1 - iphi2);
194 Int_t coldiff = TMath::Abs(ieta1 - ieta2) ;
196 // neighbours with at least common side; May 11, 2007
197 if ((coldiff==0 && TMath::Abs(rowdiff)==1) || (rowdiff==0 && TMath::Abs(coldiff)==1)) {
199 //if ((coldiff==0 && TMath::Abs(rowdiff==1)) || (rowdiff==0 && TMath::Abs(coldiff==1)) || (TMath::Abs(rowdiff)==1 && TMath::Abs(coldiff==1))) rv = 1;
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);
213 //____________________________________________________________________________
214 void AliEMCALClusterizerv1::MakeClusters()
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
218 // Mar 03, 2007 by PAI
220 if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
224 // Set up TObjArray with pointers to digits to work on
225 TObjArray *digitsC = new TObjArray();
226 TIter nextdigit(fDigitsArr);
227 AliEMCALDigit *digit;
228 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
229 digitsC->AddLast(digit);
232 double e = 0.0, ehs = 0.0;
233 TIter nextdigitC(digitsC);
234 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
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
237 digitsC->Remove(digit);
241 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n",
242 fDigitsArr->GetEntries(),fMinECut,ehs));
246 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
247 TArrayI clusterECAdigitslist(fDigitsArr->GetEntries());
249 if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()) > fECAClusteringThreshold ) ){
250 // start a new Tower RecPoint
251 if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ;
253 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
254 fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
255 recPoint = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(fNumberOfECAClusters)) ;
257 fNumberOfECAClusters++ ;
259 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
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) ;
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
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) ;
284 } // scan over digits
285 } // scan over digits already in cluster
287 AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
289 else AliFatal("Null recpoint in array!");
295 AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast()));