Initial commintment
[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)
05a92d59 19// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
20// of new IO (à la PHOS)
f1487f22 21// Mar 2007, Aleksei Pavlinov - new algoritmh of pseudo clusters
483b0559 22//////////////////////////////////////////////////////////////////////////////
23// Clusterization class. Performs clusterization (collects neighbouring active cells) and
24// unfolds the clusters having several local maxima.
25// Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
26// EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
27// parameters including input digits branch title, thresholds etc.)
28// This TTask is normally called from Reconstructioner, but can as well be used in
29// standalone mode.
30// Use Case:
31// root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")
32// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
33// //reads gAlice from header file "..."
34// root [1] cl->ExecuteTask()
35// //finds RecPoints in all events stored in galice.root
36// root [2] cl->SetDigitsBranch("digits2")
37// //sets another title for Digitis (input) branch
38// root [3] cl->SetRecPointsBranch("recp2")
39// //sets another title four output branches
40// root [4] cl->SetTowerLocalMaxCut(0.03)
41// //set clusterization parameters
42// root [5] cl->ExecuteTask("deb all time")
43// //once more finds RecPoints options are
44// // deb - print number of found rec points
45// // deb all - print number of found RecPoints and some their characteristics
46// // time - print benchmarking results
47
48// --- ROOT system ---
7dd95fbf 49#include <cassert>
483b0559 50
1d59832c 51class TROOT;
e52475ed 52#include <TH1.h>
53#include <TFile.h>
1d59832c 54class TFolder;
e52475ed 55#include <TMath.h>
56#include <TMinuit.h>
57#include <TTree.h>
1d59832c 58class TSystem;
e52475ed 59#include <TBenchmark.h>
60#include <TBrowser.h>
a1e17193 61#include <TROOT.h>
1d59832c 62
483b0559 63// --- Standard library ---
64
173558f2 65
483b0559 66// --- AliRoot header files ---
5dee926e 67#include "AliRunLoader.h"
68#include "AliRun.h"
1d59832c 69#include "AliESD.h"
5dee926e 70#include "AliEMCALLoader.h"
483b0559 71#include "AliEMCALClusterizerv1.h"
70a93198 72#include "AliEMCALRecPoint.h"
483b0559 73#include "AliEMCALDigit.h"
74#include "AliEMCALDigitizer.h"
483b0559 75#include "AliEMCAL.h"
05a92d59 76#include "AliEMCALGeometry.h"
98e9578e 77#include "AliEMCALRawUtils.h"
e52475ed 78#include "AliEMCALHistoUtilities.h"
1bd98442 79#include "AliCDBManager.h"
1d59832c 80
81class AliCDBStorage;
1bd98442 82#include "AliCDBEntry.h"
483b0559 83
84ClassImp(AliEMCALClusterizerv1)
1963b290 85
483b0559 86//____________________________________________________________________________
18a21c7c 87AliEMCALClusterizerv1::AliEMCALClusterizerv1()
88 : AliEMCALClusterizer(),
89 fHists(0),fPointE(0),fPointL1(0),fPointL2(0),
90 fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0),
91 fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0),
92 fDefaultInit(kTRUE),
93 fToUnfold(kFALSE),
94 fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0),
95 fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
98e9578e 96 fECAW0(0.),fRecPointsInRun(0),fTimeCut(0.),fMinECut(0.)
483b0559 97{
98 // default ctor (to be used mainly by Streamer)
99
839828a6 100 InitParameters() ;
54b24bce 101 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
a5c60732 102 fGeom->GetTransformationForSM(); // Global <-> Local
483b0559 103}
104
105//____________________________________________________________________________
88cb7938 106AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
18a21c7c 107 : AliEMCALClusterizer(alirunFileName, eventFolderName),
108 fHists(0),fPointE(0),fPointL1(0),fPointL2(0),
109 fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0),
110 fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0),
111 fDefaultInit(kFALSE),
112 fToUnfold(kFALSE),
113 fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0),
114 fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.),
98e9578e 115 fECAW0(0.),fRecPointsInRun(0),fTimeCut(0.),fMinECut(0.)
483b0559 116{
117 // ctor with the indication of the file where header Tree and digits Tree are stored
118
839828a6 119 InitParameters() ;
483b0559 120 Init() ;
483b0559 121}
05a92d59 122
483b0559 123//____________________________________________________________________________
18a21c7c 124AliEMCALClusterizerv1::AliEMCALClusterizerv1(const AliEMCALClusterizerv1& clus)
125 : AliEMCALClusterizer(clus),
126 fHists(clus.fHists),
127 fPointE(clus.fPointE),
128 fPointL1(clus.fPointL1),
129 fPointL2(clus.fPointL2),
130 fPointDis(clus.fPointDis),
131 fPointMult(clus.fPointMult),
132 fDigitAmp(clus.fDigitAmp),
133 fMaxE(clus.fMaxE),
134 fMaxL1(clus.fMaxL1),
135 fMaxL2(clus.fMaxL2),
136 fMaxDis(clus.fMaxDis),
137 fGeom(clus.fGeom),
138 fDefaultInit(clus.fDefaultInit),
139 fToUnfold(clus.fToUnfold),
140 fNumberOfECAClusters(clus.fNumberOfECAClusters),
141 fNTowerInGroup(clus.fNTowerInGroup),
142 fCalibData(clus.fCalibData),
143 fADCchannelECA(clus.fADCchannelECA),
144 fADCpedestalECA(clus.fADCpedestalECA),
145 fECAClusteringThreshold(clus.fECAClusteringThreshold),
146 fECALocMaxCut(clus.fECALocMaxCut),
147 fECAW0(clus.fECAW0),
148 fRecPointsInRun(clus.fRecPointsInRun),
98e9578e 149 fTimeCut(clus.fTimeCut),
18a21c7c 150 fMinECut(clus.fMinECut)
0a4cb131 151{
152 //copy ctor
0a4cb131 153}
154
155//____________________________________________________________________________
483b0559 156 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
157{
ef305168 158 // dtor
ef305168 159}
160
161//____________________________________________________________________________
162const TString AliEMCALClusterizerv1::BranchName() const
9859bfc0 163{
88cb7938 164 return GetName();
483b0559 165}
839828a6 166
483b0559 167//____________________________________________________________________________
1bd98442 168Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId)
9859bfc0 169{
1bd98442 170
171 // Convert digitized amplitude into energy.
172 // Calibration parameters are taken from calibration data base for raw data,
173 // or from digitizer parameters for simulated data.
174
175 if(fCalibData){
1bd98442 176
a5c60732 177 if (fGeom==0)
1bd98442 178 AliFatal("Did not get geometry from EMCALLoader") ;
179
180 Int_t iSupMod = -1;
2bb3725c 181 Int_t nModule = -1;
1bd98442 182 Int_t nIphi = -1;
183 Int_t nIeta = -1;
184 Int_t iphi = -1;
185 Int_t ieta = -1;
186
2bb3725c 187 Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ;
1d46d1f6 188 if(!bCell) {
189 fGeom->PrintGeometry();
190 Error("Calibrate()"," Wrong cell id number : %i", AbsId);
191 assert(0);
192 }
54b24bce 193
2bb3725c 194 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
1bd98442 195
196 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
197 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
54b24bce 198
199 return -fADCpedestalECA + amp * fADCchannelECA ;
1bd98442 200
201 }
202 else //Return energy with default parameters if calibration is not available
203 return -fADCpedestalECA + amp * fADCchannelECA ;
204
483b0559 205}
05a92d59 206
483b0559 207//____________________________________________________________________________
208void AliEMCALClusterizerv1::Exec(Option_t * option)
209{
98e9578e 210 // Steering method to perform clusterization for the current event
211 // in AliEMCALLoader
483b0559 212
483b0559 213 if(strstr(option,"tim"))
214 gBenchmark->Start("EMCALClusterizer");
215
216 if(strstr(option,"print"))
217 Print("") ;
98e9578e 218
5dee926e 219 AliRunLoader *rl = AliRunLoader::GetRunLoader();
220 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(rl->GetDetectorLoader("EMCAL"));
88cb7938 221
1bd98442 222 //Get calibration parameters from file or digitizer default values.
223 GetCalibrationParameters() ;
224
483b0559 225
98e9578e 226 fNumberOfECAClusters = 0;
483b0559 227
98e9578e 228 if(strstr(option,"pseudo"))
229 MakeClusters("pseudo") ; //both types
230 else
231 MakeClusters("") ; //only the real clusters
a5c60732 232
98e9578e 233 if(fToUnfold)
234 MakeUnfolding() ;
ed611565 235
98e9578e 236 WriteRecPoints() ;
483b0559 237
98e9578e 238 if(strstr(option,"deb") || strstr(option,"all"))
239 PrintRecPoints(option) ;
483b0559 240
98e9578e 241 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",emcalLoader->RecPoints()->GetEntriesFast()));
483b0559 242
98e9578e 243 //increment the total number of recpoints per run
244 fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ;
88cb7938 245
483b0559 246 if(strstr(option,"tim")){
247 gBenchmark->Stop("EMCALClusterizer");
98e9578e 248 printf("Exec took %f seconds for Clusterizing",
249 gBenchmark->GetCpuTime("EMCALClusterizer"));
250 }
483b0559 251}
252
253//____________________________________________________________________________
70a93198 254Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
483b0559 255 Int_t nPar, Float_t * fitparameters) const
256{
257 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
258 // The initial values for fitting procedure are set equal to the positions of local maxima.
259 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
260
5dee926e 261 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
262 TClonesArray *digits = emcalLoader->Digits();
5dee926e 263
483b0559 264 gMinuit->mncler(); // Reset Minuit's list of paramters
265 gMinuit->SetPrintLevel(-1) ; // No Printout
266 gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
267 // To set the address of the minimization function
483b0559 268 TList * toMinuit = new TList();
269 toMinuit->AddAt(emcRP,0) ;
270 toMinuit->AddAt(digits,1) ;
271
272 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
273
274 // filling initial values for fit parameters
275 AliEMCALDigit * digit ;
276
277 Int_t ierflg = 0;
278 Int_t index = 0 ;
279 Int_t nDigits = (Int_t) nPar / 3 ;
280
281 Int_t iDigit ;
282
483b0559 283 for(iDigit = 0; iDigit < nDigits; iDigit++){
a0636361 284 digit = maxAt[iDigit];
483b0559 285
483b0559 286 Float_t x = 0.;
287 Float_t z = 0.;
dc7da436 288 // have to be tune for TRD1; May 31,06
289 // Int_t relid[2] ;
290 // fGeom->AbsToRelNumbering(digit->GetId(), relid) ; // obsolete method
291 // fGeom->PosInAlice(relid, x, z) ; // obsolete method
483b0559 292
293 Float_t energy = maxAtEnergy[iDigit] ;
294
295 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
296 index++ ;
297 if(ierflg != 0){
9859bfc0 298 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ;
483b0559 299 return kFALSE;
300 }
301 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
302 index++ ;
303 if(ierflg != 0){
9859bfc0 304 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ;
483b0559 305 return kFALSE;
306 }
307 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
308 index++ ;
309 if(ierflg != 0){
9859bfc0 310 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ;
483b0559 311 return kFALSE;
312 }
313 }
314
315 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
316 // depends on it.
317 Double_t p1 = 1.0 ;
318 Double_t p2 = 0.0 ;
319
320 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
321 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
322 gMinuit->SetMaxIterations(5);
323 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
324
325 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
326
327 if(ierflg == 4){ // Minimum not found
9859bfc0 328 Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ;
483b0559 329 return kFALSE ;
330 }
331 for(index = 0; index < nPar; index++){
332 Double_t err ;
333 Double_t val ;
334 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
335 fitparameters[index] = val ;
336 }
337
338 delete toMinuit ;
339 return kTRUE;
340
341}
342
343//____________________________________________________________________________
344void AliEMCALClusterizerv1::GetCalibrationParameters()
345{
1bd98442 346 // Set calibration parameters:
347 // if calibration database exists, they are read from database,
348 // otherwise, they are taken from digitizer.
349 //
350 // It is a user responsilibity to open CDB before reconstruction,
351 // for example:
352 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
353
354 //Check if calibration is stored in data base
54b24bce 355
356 AliEMCALLoader *emcalLoader =
357 dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
358
359 fCalibData =emcalLoader->CalibData();
360
1bd98442 361 if(!fCalibData)
362 {
363 //If calibration is not available use default parameters
364 //Loader
1bd98442 365 if ( !emcalLoader->Digitizer() )
366 emcalLoader->LoadDigitizer();
367 AliEMCALDigitizer * dig = dynamic_cast<AliEMCALDigitizer*>(emcalLoader->Digitizer());
368
369 fADCchannelECA = dig->GetECAchannel() ;
370 fADCpedestalECA = dig->GetECApedestal();
371 }
483b0559 372}
05a92d59 373
483b0559 374//____________________________________________________________________________
375void AliEMCALClusterizerv1::Init()
376{
377 // Make all memory allocations which can not be done in default constructor.
378 // Attach the Clusterizer task to the list of EMCAL tasks
379
5dee926e 380 AliRunLoader *rl = AliRunLoader::GetRunLoader();
98e9578e 381 if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL"))
382 fGeom = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry();
383 else
384 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName());
385
a5c60732 386 fGeom->GetTransformationForSM(); // Global <-> Local
387 AliInfo(Form("geom 0x%x",fGeom));
05a92d59 388
483b0559 389 if(!gMinuit)
390 gMinuit = new TMinuit(100) ;
a5c60732 391
392 fHists = BookHists();
483b0559 393}
394
395//____________________________________________________________________________
839828a6 396void AliEMCALClusterizerv1::InitParameters()
fdebddeb 397{
398 // Initializes the parameters for the Clusterizer
399 fNumberOfECAClusters = 0;
1963b290 400
a5c60732 401 fNTowerInGroup = 36; //Produces maximum of 80 pseudoclusters per event
402
990491dc 403 fECAClusteringThreshold = 0.5; // Best value for 2 GeV gamma merged with Ideal HIJING. Retune later?
e52475ed 404 fECALocMaxCut = 0.03; // ??
1963b290 405
1d46d1f6 406 fECAW0 = 4.5;
98e9578e 407 fTimeCut = 300e-9 ; // 300 ns time cut (to be tuned)
839828a6 408 fToUnfold = kFALSE ;
fdebddeb 409 fRecPointsInRun = 0 ;
990491dc 410 fMinECut = 0.45; // Best value for 2 GeV gamma merged with Ideal HIJING. Retune later?
1bd98442 411
412 fCalibData = 0 ;
3a8be91c 413
414 // If reconstruction parameters are found in OCDB, take them from it
415
416 AliEMCALLoader *emcalLoader =
417 dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
418 AliEMCALRecParam *recParamDB = emcalLoader->RecParam();
419 if (recParamDB != 0) {
420 fECAClusteringThreshold = recParamDB->GetClusteringThreshold();
421 fECAW0 = recParamDB->GetW0();
422 fMinECut = recParamDB->GetMinECut();
423 AliDebug(1,Form("Reconstruction parameters were taken from OCDB: fECAClusteringThreshold=%.3f, fECAW=%.3f, fMinECut=%.3f",
424 fECAClusteringThreshold,fECAW0,fMinECut));
425 }
426
839828a6 427}
428
429//____________________________________________________________________________
e52475ed 430Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const
483b0559 431{
1d46d1f6 432 // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching
483b0559 433 // = 1 are neighbour
1d46d1f6 434 // = 2 is in different SM; continue searching
483b0559 435 // neighbours are defined as digits having at least a common vertex
436 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
437 // which is compared to a digit (d2) not yet in a cluster
438
e52475ed 439 static Int_t rv;
2bb3725c 440 static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
441 static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
e52475ed 442 static Int_t rowdiff, coldiff;
443 rv = 0 ;
483b0559 444
2bb3725c 445 fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
446 fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
e52475ed 447 if(nSupMod1 != nSupMod2) return 2; // different SM
448
2bb3725c 449 fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1);
450 fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2);
e52475ed 451
452 rowdiff = TMath::Abs(iphi1 - iphi2);
453 coldiff = TMath::Abs(ieta1 - ieta2) ;
70a93198 454
3d99d6d5 455 // neighbours with at least commom side; May 11, 2007
456 if ((coldiff==0 && abs(rowdiff)==1) || (rowdiff==0 && abs(coldiff)==1)) rv = 1;
70a93198 457
e52475ed 458 if (gDebug == 2 && rv==1)
459 printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
460 rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
ed611565 461
483b0559 462 return rv ;
463}
464
88cb7938 465//____________________________________________________________________________
a5c60732 466Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) const
467{
468 // Tells whether two digits fall within the same supermodule and
469 // tower grouping. The number of towers in a group is controlled by
470 // the parameter nTowersInGroup
471 // = 0 are not in same group but continue searching
472 // = 1 same group
473 // = 2 is in different SM, quit from searching
474 // = 3 different tower group, quit from searching
475 //
476 // The order of d1 and d2 is important: first (d1) should be a digit
477 // already in a cluster which is compared to a digit (d2) not yet in a cluster
478
479 //JLK Question: does the quit from searching assume that the digits
480 //are ordered, so that once you are in a different SM, you'll not
481 //find another in the list that will match? How about my TowerGroup search?
482
483 static Int_t rv;
2bb3725c 484 static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0;
485 static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0;
a5c60732 486 static Int_t towerGroup1 = -1, towerGroup2 = -1;
487 rv = 0 ;
488
2bb3725c 489 fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1);
490 fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2);
a5c60732 491 if(nSupMod1 != nSupMod2) return 2; // different SM
492
2bb3725c 493 static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInModule();
a5c60732 494
495 //figure out which tower grouping each digit belongs to
496 for(int it = 0; it < nTowerInSM/fNTowerInGroup; it++) {
2bb3725c 497 if(nModule1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it;
498 if(nModule2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it;
a5c60732 499 }
500 if(towerGroup1 != towerGroup2) return 3; //different Towergroup
501
502 //same SM, same towergroup, we're happy
503 if(towerGroup1 == towerGroup2 && towerGroup2 >= 0)
504 rv = 1;
505
506 if (gDebug == 2 && rv==1)
507 printf("AreInGroup: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n",
508 rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2);
509
510 return rv ;
511}
483b0559 512
483b0559 513//____________________________________________________________________________
9e5d2067 514void AliEMCALClusterizerv1::WriteRecPoints()
483b0559 515{
516
517 // Creates new branches with given title
518 // fills and writes into TreeR.
519
5dee926e 520 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
ed611565 521
5dee926e 522 TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
ed611565 523
5dee926e 524 TClonesArray * digits = emcalLoader->Digits() ;
525 TTree * treeR = emcalLoader->TreeR();
526 if ( treeR==0 ) {
527 emcalLoader->MakeRecPointsContainer();
528 treeR = emcalLoader->TreeR();
529 }
98e9578e 530 else if (treeR->GetEntries() > 0) {
531 Warning("WriteRecPoints","RecPoints already exist in output file. New Recpoitns will not be visible.");
532 }
05a92d59 533 Int_t index ;
ed611565 534
ed611565 535 //Evaluate position, dispersion and other RecPoint properties for EC section
98e9578e 536 for(index = 0; index < aECARecPoints->GetEntries(); index++) {
537 if (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index))->GetClusterType() != AliESDCaloCluster::kPseudoCluster)
538 dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index))->EvalAll(fECAW0,digits) ;
539 }
ed611565 540
88cb7938 541 aECARecPoints->Sort() ;
483b0559 542
e52475ed 543 for(index = 0; index < aECARecPoints->GetEntries(); index++) {
70a93198 544 (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
e52475ed 545 (dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(index)))->Print();
546 }
483b0559 547
483b0559 548 Int_t bufferSize = 32000 ;
9859bfc0 549 Int_t splitlevel = 0 ;
ed611565 550
551 //EC section branch
5dee926e 552
553 TBranch * branchECA = 0;
554 if ((branchECA = treeR->GetBranch("EMCALECARP")))
555 branchECA->SetAddress(&aECARecPoints);
556 else
557 treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
558 treeR->Fill() ;
483b0559 559
5dee926e 560 emcalLoader->WriteRecPoints("OVERWRITE");
483b0559 561}
562
563//____________________________________________________________________________
a5c60732 564void AliEMCALClusterizerv1::MakeClusters(char* option)
483b0559 565{
566 // Steering method to construct the clusters stored in a list of Reconstructed Points
567 // A cluster is defined as a list of neighbour digits
f1487f22 568 // Mar 03, 2007 by PAI
ed611565 569
f1487f22 570 if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader");
ed611565 571
f1487f22 572 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
5dee926e 573 TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
5dee926e 574 aECARecPoints->Clear();
ed611565 575
f1487f22 576 TClonesArray *digits = emcalLoader->Digits();
98e9578e 577
578 // Set up TObjArray with pointers to digits to work on
579 TObjArray *digitsC = new TObjArray();
580 TIter nextdigit(digits);
581 AliEMCALDigit *digit;
582 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextdigit())) ) {
583 digitsC->AddLast(digit);
584 }
a5c60732 585
f1487f22 586 //Start with pseudoclusters, if option
587 if(strstr(option,"pseudo")) {
588 //
589 // New algorithm : will be created one pseudo cluster per module
590 //
591 AliDebug(1,Form("Pseudo clustering #digits : %i ",digits->GetEntries()));
592
593 AliEMCALRecPoint *recPoints[12]; // max size is 12 : see fGeom->GetNumberOfSuperModules();
594 for(int i=0; i<12; i++) recPoints[i] = 0;
98e9578e 595 TIter nextdigitC(digitsC) ;
a5c60732 596
98e9578e 597 // PseudoClusterization starts
f1487f22 598 int nSM = 0; // # of SM
98e9578e 599 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
a5c60732 600 if(fGeom->CheckAbsCellId(digit->GetId()) ) { //Is this an EMCAL digit? Just maing sure...
f1487f22 601 nSM = fGeom->GetSuperModuleNumber(digit->GetId());
602 if(recPoints[nSM] == 0) {
603 recPoints[nSM] = new AliEMCALRecPoint(Form("PC%2.2i", nSM));
604 recPoints[nSM]->SetClusterType(AliESDCaloCluster::kPseudoCluster);
605 }
606 recPoints[nSM]->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId()));
a5c60732 607 }
608 }
98e9578e 609 fNumberOfECAClusters = 0;
f1487f22 610 for(int i=0; i<fGeom->GetNumberOfSuperModules(); i++) { // put non empty rec.points to container
611 if(recPoints[i]) aECARecPoints->AddAt(recPoints[i], fNumberOfECAClusters++);
612 }
613 AliDebug(1,Form(" Number of PC %d ", fNumberOfECAClusters));
a5c60732 614 }
615
f1487f22 616 //
617 // Now do real clusters
618 //
a5c60732 619
1d46d1f6 620 double e = 0.0, ehs = 0.0;
98e9578e 621 TIter nextdigitC(digitsC);
622
1d46d1f6 623 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // clean up digits
1bd98442 624 e = Calibrate(digit->GetAmp(), digit->GetId());
e52475ed 625 AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp());
626 AliEMCALHistoUtilities::FillH1(fHists, 11, e);
98e9578e 627 if ( e < fMinECut || digit->GetTimeR() > fTimeCut )
628 digitsC->Remove(digit);
629 else
630 ehs += e;
1d59832c 631 }
632 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %d\n",
633 digits->GetEntries(),fMinECut,ehs));
98e9578e 634
1d46d1f6 635 nextdigitC.Reset();
98e9578e 636
1d46d1f6 637 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigitC())) ) { // scan over the list of digitsC
f1487f22 638 TArrayI clusterECAdigitslist(digits->GetEntries());
e52475ed 639
a5c60732 640 if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold ) ){
98e9578e 641 // start a new Tower RecPoint
e52475ed 642 if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
98e9578e 643 AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ;
644 aECARecPoints->AddAt(recPoint, fNumberOfECAClusters) ;
e52475ed 645 recPoint = dynamic_cast<AliEMCALRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
4aa81449 646 fNumberOfECAClusters++ ;
a5c60732 647
1d59832c 648 recPoint->SetClusterType(AliESDCaloCluster::kClusterv1);
a5c60732 649
1bd98442 650 recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ;
98e9578e 651 TObjArray clusterDigits;
652 clusterDigits.AddLast(digit);
4aa81449 653 digitsC->Remove(digit) ;
f1487f22 654
655 AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(),
656 Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold));
483b0559 657
98e9578e 658 // Grow cluster by finding neighbours
659 TIter nextClusterDigit(&clusterDigits);
660 while ( (digit = dynamic_cast<AliEMCALDigit*>(nextClusterDigit())) ) { // scan over digits in cluster
661 TIter nextdigitN(digitsC);
662 AliEMCALDigit *digitN = 0; // digi neighbor
663 while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours
664 if (AreNeighbours(digit, digitN)==1) { // call (digit,digitN) in THAT oder !!!!!
1d46d1f6 665 recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ;
98e9578e 666 clusterDigits.AddLast(digitN) ;
1d46d1f6 667 digitsC->Remove(digitN) ;
1d46d1f6 668 } // if(ineb==1)
98e9578e 669 } // scan over digits
e52475ed 670 } // scan over digits already in cluster
98e9578e 671 if(recPoint)
672 AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy()));
673 } // If seed found
1d59832c 674 } // while digit
98e9578e 675
483b0559 676 delete digitsC ;
a5c60732 677
5dee926e 678 AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast()));
483b0559 679}
680
fdebddeb 681void AliEMCALClusterizerv1::MakeUnfolding() const
483b0559 682{
683 Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
483b0559 684}
685
686//____________________________________________________________________________
687Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r)
688{
689 // Shape of the shower (see EMCAL TDR)
690 // If you change this function, change also the gradient evaluation in ChiSquare()
691
692 Double_t r4 = r*r*r*r ;
693 Double_t r295 = TMath::Power(r, 2.95) ;
694 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
695 return shape ;
696}
697
698//____________________________________________________________________________
70a93198 699void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/,
9e5d2067 700 Int_t /*nMax*/,
701 AliEMCALDigit ** /*maxAt*/,
fdebddeb 702 Float_t * /*maxAtEnergy*/) const
483b0559 703{
704 // Performs the unfolding of a cluster with nMax overlapping showers
483b0559 705
9859bfc0 706 Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
173558f2 707
483b0559 708}
709
710//_____________________________________________________________________________
9e5d2067 711void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/,
712 Double_t & /*fret*/,
713 Double_t * /*x*/, Int_t /*iflag*/)
483b0559 714{
715 // Calculates the Chi square for the cluster unfolding minimization
716 // Number of parameters, Gradient, Chi squared, parameters, what to do
483b0559 717
9859bfc0 718 ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
483b0559 719}
483b0559 720//____________________________________________________________________________
9e5d2067 721void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const
483b0559 722{
723 // Print clusterizer parameters
724
9859bfc0 725 TString message("\n") ;
726
483b0559 727 if( strcmp(GetName(), "") !=0 ){
728
729 // Print parameters
730
fdebddeb 731 TString taskName(GetName()) ;
483b0559 732 taskName.ReplaceAll(Version(), "") ;
9859bfc0 733
fdebddeb 734 printf("--------------- ");
735 printf(taskName.Data()) ;
736 printf(" ");
737 printf(GetTitle()) ;
738 printf("-----------\n");
739 printf("Clusterizing digits from the file: ");
740 printf(taskName.Data());
741 printf("\n Branch: ");
742 printf(GetName());
743 printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
b481a360 744 printf("\n ECA Logarithmic weight = %f", fECAW0);
483b0559 745 if(fToUnfold)
fdebddeb 746 printf("\nUnfolding on\n");
483b0559 747 else
fdebddeb 748 printf("\nUnfolding off\n");
483b0559 749
fdebddeb 750 printf("------------------------------------------------------------------");
483b0559 751 }
752 else
fdebddeb 753 printf("AliEMCALClusterizerv1 not initialized ") ;
483b0559 754}
173558f2 755
483b0559 756//____________________________________________________________________________
757void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
758{
b481a360 759 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
5dee926e 760 AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL"));
761 TObjArray * aECARecPoints = emcalLoader->RecPoints() ;
762
1d46d1f6 763 if(strstr(option,"deb")) {
764 printf("PrintRecPoints: Clusterization result:") ;
ed611565 765
1d46d1f6 766 printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ;
767 printf(" Found %d ECA Rec Points\n ",
fdebddeb 768 aECARecPoints->GetEntriesFast()) ;
1d46d1f6 769 }
483b0559 770
88cb7938 771 fRecPointsInRun += aECARecPoints->GetEntriesFast() ;
11f9c5ff 772
483b0559 773 if(strstr(option,"all")) {
1d46d1f6 774 if(strstr(option,"deb")) {
775 printf("\n-----------------------------------------------------------------------\n") ;
776 printf("Clusters in ECAL section\n") ;
777 printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
778 }
779 Int_t index =0;
1963b290 780 Float_t maxE=0;
781 Float_t maxL1=0;
782 Float_t maxL2=0;
783 Float_t maxDis=0;
784
1d46d1f6 785 AliEMCALHistoUtilities::FillH1(fHists, 12, double(aECARecPoints->GetEntries()));
786
88cb7938 787 for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
70a93198 788 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(aECARecPoints->At(index)) ;
483b0559 789 TVector3 globalpos;
1963b290 790 //rp->GetGlobalPosition(globalpos);
ed611565 791 TVector3 localpos;
792 rp->GetLocalPosition(localpos);
483b0559 793 Float_t lambda[2];
794 rp->GetElipsAxis(lambda);
795 Int_t * primaries;
796 Int_t nprimaries;
797 primaries = rp->GetPrimaries(nprimaries);
1d46d1f6 798 if(strstr(option,"deb"))
70a93198 799 printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
800 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
ed611565 801 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
802 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
1963b290 803 /////////////
804 if(rp->GetEnergy()>maxE){
805 maxE=rp->GetEnergy();
806 maxL1=lambda[0];
807 maxL2=lambda[1];
808 maxDis=rp->GetDispersion();
809 }
a5c60732 810 fPointE->Fill(rp->GetEnergy());
811 fPointL1->Fill(lambda[0]);
812 fPointL2->Fill(lambda[1]);
813 fPointDis->Fill(rp->GetDispersion());
814 fPointMult->Fill(rp->GetMultiplicity());
1963b290 815 /////////////
1d46d1f6 816 if(strstr(option,"deb")){
817 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
818 printf("%d ", primaries[iprimary] ) ;
819 }
820 }
483b0559 821 }
1963b290 822
a5c60732 823 fMaxE->Fill(maxE);
824 fMaxL1->Fill(maxL1);
825 fMaxL2->Fill(maxL2);
826 fMaxDis->Fill(maxDis);
1963b290 827
1d46d1f6 828 if(strstr(option,"deb"))
ed611565 829 printf("\n-----------------------------------------------------------------------\n");
483b0559 830 }
831}
e52475ed 832TList* AliEMCALClusterizerv1::BookHists()
1963b290 833{
14ce0a6e 834 //set up histograms for monitoring clusterizer performance
835
e52475ed 836 gROOT->cd();
837
1d46d1f6 838 fPointE = new TH1F("00_pointE","point energy", 2000, 0.0, 150.);
839 fPointL1 = new TH1F("01_pointL1","point L1", 1000, 0.0, 3.);
840 fPointL2 = new TH1F("02_pointL2","point L2", 1000, 0.0, 3.);
841 fPointDis = new TH1F("03_pointDisp","point dispersion", 1000, 0.0, 10.);
842 fPointMult = new TH1F("04_pointMult","#cell in point(cluster)", 101, -0.5, 100.5);
843 fDigitAmp = new TH1F("05_digitAmp","Digit Amplitude", 2000, 0.0, 5000.);
844 fMaxE = new TH1F("06_maxE","Max point energy", 2000, 0.0, 150.);
845 fMaxL1 = new TH1F("07_maxL1","Largest (first) of eigenvalue of covariance matrix", 1000, 0.0, 3.);
846 fMaxL2 = new TH1F("08_maxL2","Smalest (second) of eigenvalue of covariace matrix", 1000, 0.0, 3.);
847 fMaxDis = new TH1F("09_maxDis","Point dispersion", 1000, 0.0, 10.); // 9
e52475ed 848 //
1d46d1f6 849 new TH1F("10_adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5); // 10
850 new TH1F("11_energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.); // 11
851 new TH1F("12_numberOfPoints","number of points(clusters)", 101, -0.5, 100.5); // 12
e52475ed 852
853 return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE);
1963b290 854}
e52475ed 855
856void AliEMCALClusterizerv1::SaveHists(const char *fn)
1963b290 857{
e52475ed 858 AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE);
1963b290 859}
860
1d46d1f6 861void AliEMCALClusterizerv1::PrintRecoInfo()
862{
863 printf(" AliEMCALClusterizerv1::PrintRecoInfo() : version %s \n", Version() );
864 TH1F *h = (TH1F*)fHists->At(12);
865 if(h) {
866 printf(" ## Multiplicity of RecPoints ## \n");
867 for(int i=1; i<=h->GetNbinsX(); i++) {
868 int nbin = int((*h)[i]);
869 int mult = int(h->GetBinCenter(i));
870 if(nbin > 0) printf(" %i : %5.5i %6.3f %% \n", mult, nbin, 100.*nbin/h->GetEntries());
871 }
872 }
873}
874
875void AliEMCALClusterizerv1::DrawLambdasHists()
876{
877 if(fMaxL1) {
878 fMaxL1->Draw();
879 if(fMaxL2) fMaxL2->Draw("same");
880 if(fMaxDis) {
881 fMaxDis->Draw("same");
882 }
883 }
884}
885
e52475ed 886void AliEMCALClusterizerv1::Browse(TBrowser* b)
1963b290 887{
e52475ed 888 if(fHists) b->Add(fHists);
1d46d1f6 889 if(fGeom) b->Add(fGeom);
e52475ed 890 TTask::Browse(b);
1963b290 891}