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