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