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