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