Transition to NewIO
[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$ */
483b0559 17/* $Log:
18 1 October 2000. Yuri Kharlov:
19 AreNeighbours()
20 PPSD upper layer is considered if number of layers>1
483b0559 21 18 October 2000. Yuri Kharlov:
22 AliEMCALClusterizerv1()
23 CPV clusterizing parameters added
483b0559 24 MakeClusters()
25 After first PPSD digit remove EMC digits only once
26*/
27//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
05a92d59 28// August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction
29// of new IO (à la PHOS)
483b0559 30//////////////////////////////////////////////////////////////////////////////
31// Clusterization class. Performs clusterization (collects neighbouring active cells) and
32// unfolds the clusters having several local maxima.
33// Results are stored in TreeR#, branches EMCALTowerRP (EMC recPoints),
34// EMCALPreShoRP (CPV RecPoints) and AliEMCALClusterizer (Clusterizer with all
35// parameters including input digits branch title, thresholds etc.)
36// This TTask is normally called from Reconstructioner, but can as well be used in
37// standalone mode.
38// Use Case:
39// root [0] AliEMCALClusterizerv1 * cl = new AliEMCALClusterizerv1("galice.root")
40// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
41// //reads gAlice from header file "..."
42// root [1] cl->ExecuteTask()
43// //finds RecPoints in all events stored in galice.root
44// root [2] cl->SetDigitsBranch("digits2")
45// //sets another title for Digitis (input) branch
46// root [3] cl->SetRecPointsBranch("recp2")
47// //sets another title four output branches
48// root [4] cl->SetTowerLocalMaxCut(0.03)
49// //set clusterization parameters
50// root [5] cl->ExecuteTask("deb all time")
51// //once more finds RecPoints options are
52// // deb - print number of found rec points
53// // deb all - print number of found RecPoints and some their characteristics
54// // time - print benchmarking results
55
56// --- ROOT system ---
57
58#include "TROOT.h"
59#include "TFile.h"
60#include "TFolder.h"
61#include "TMath.h"
62#include "TMinuit.h"
63#include "TTree.h"
64#include "TSystem.h"
65#include "TBenchmark.h"
66
67// --- Standard library ---
68
173558f2 69
483b0559 70// --- AliRoot header files ---
88cb7938 71#include "AliEMCALGetter.h"
483b0559 72#include "AliEMCALClusterizerv1.h"
88cb7938 73#include "AliEMCALTowerRecPoint.h"
483b0559 74#include "AliEMCALDigit.h"
75#include "AliEMCALDigitizer.h"
483b0559 76#include "AliEMCAL.h"
05a92d59 77#include "AliEMCALGeometry.h"
483b0559 78
79ClassImp(AliEMCALClusterizerv1)
80
81//____________________________________________________________________________
82 AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer()
83{
84 // default ctor (to be used mainly by Streamer)
85
839828a6 86 InitParameters() ;
ef305168 87 fDefaultInit = kTRUE ;
483b0559 88}
89
90//____________________________________________________________________________
88cb7938 91AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName)
92:AliEMCALClusterizer(alirunFileName, eventFolderName)
483b0559 93{
94 // ctor with the indication of the file where header Tree and digits Tree are stored
95
839828a6 96 InitParameters() ;
483b0559 97 Init() ;
05a92d59 98 fDefaultInit = kFALSE ;
483b0559 99
100}
05a92d59 101
483b0559 102//____________________________________________________________________________
103 AliEMCALClusterizerv1::~AliEMCALClusterizerv1()
104{
ef305168 105 // dtor
839828a6 106
ef305168 107}
108
109//____________________________________________________________________________
110const TString AliEMCALClusterizerv1::BranchName() const
9859bfc0 111{
88cb7938 112 return GetName();
113
483b0559 114}
839828a6 115
483b0559 116//____________________________________________________________________________
ed611565 117Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t where) const
9859bfc0 118{
119 //To be replased later by the method, reading individual parameters from the database
ed611565 120 // where = 0 == PRE ; where = 1 == ECAL ; where = 2 == HCAL
121 if ( where == 0 ) // calibrate as PRE section
122 return -fADCpedestalPRE + amp * fADCchannelPRE ;
88cb7938 123 else if (where == 1) //calibrate as ECA section
124 return -fADCpedestalECA + amp * fADCchannelECA ;
125 else if (where == 2) //calibrate as HCA section
126 return -fADCpedestalHCA + amp * fADCchannelHCA ;
ed611565 127 else
128 Fatal("Calibrate", "Something went wrong!") ;
129 return -9999999. ;
483b0559 130}
05a92d59 131
483b0559 132//____________________________________________________________________________
133void AliEMCALClusterizerv1::Exec(Option_t * option)
134{
135 // Steering method
136
483b0559 137 if(strstr(option,"tim"))
138 gBenchmark->Start("EMCALClusterizer");
139
140 if(strstr(option,"print"))
141 Print("") ;
142
88cb7938 143 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
144
05a92d59 145 Int_t nevents = gime->MaxEvent() ;
483b0559 146 Int_t ievent ;
147
148 for(ievent = 0; ievent < nevents; ievent++){
149
05a92d59 150 gime->Event(ievent,"D") ;
151
88cb7938 152 GetCalibrationParameters() ;
483b0559 153
88cb7938 154 fNumberOfPREClusters = fNumberOfECAClusters = fNumberOfHCAClusters = 0 ;
05a92d59 155
483b0559 156 MakeClusters() ;
ed611565 157
483b0559 158 if(fToUnfold)
159 MakeUnfolding() ;
160
161 WriteRecPoints(ievent) ;
162
163 if(strstr(option,"deb"))
164 PrintRecPoints(option) ;
165
88cb7938 166 //increment the total number of recpoints per run
ed611565 167 fRecPointsInRun += gime->PRERecPoints()->GetEntriesFast() ;
88cb7938 168 fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ;
169 fRecPointsInRun += gime->HCARecPoints()->GetEntriesFast() ;
170 }
483b0559 171
88cb7938 172 Unload();
173
483b0559 174 if(strstr(option,"tim")){
175 gBenchmark->Stop("EMCALClusterizer");
9859bfc0 176 Info("Exec", "took %f seconds for Clusterizing %f seconds per event",
177 gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nevents ) ;
483b0559 178 }
483b0559 179}
180
181//____________________________________________________________________________
a0636361 182Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy,
483b0559 183 Int_t nPar, Float_t * fitparameters) const
184{
185 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
186 // The initial values for fitting procedure are set equal to the positions of local maxima.
187 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
188
88cb7938 189 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
483b0559 190 TClonesArray * digits = gime->Digits() ;
191
192
193 gMinuit->mncler(); // Reset Minuit's list of paramters
194 gMinuit->SetPrintLevel(-1) ; // No Printout
195 gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ;
196 // To set the address of the minimization function
483b0559 197 TList * toMinuit = new TList();
198 toMinuit->AddAt(emcRP,0) ;
199 toMinuit->AddAt(digits,1) ;
200
201 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
202
203 // filling initial values for fit parameters
204 AliEMCALDigit * digit ;
205
206 Int_t ierflg = 0;
207 Int_t index = 0 ;
208 Int_t nDigits = (Int_t) nPar / 3 ;
209
210 Int_t iDigit ;
211
212 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
213
214 for(iDigit = 0; iDigit < nDigits; iDigit++){
a0636361 215 digit = maxAt[iDigit];
483b0559 216
217 Int_t relid[4] ;
218 Float_t x = 0.;
219 Float_t z = 0.;
220 geom->AbsToRelNumbering(digit->GetId(), relid) ;
221 geom->PosInAlice(relid, x, z) ;
222
223 Float_t energy = maxAtEnergy[iDigit] ;
224
225 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
226 index++ ;
227 if(ierflg != 0){
9859bfc0 228 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ;
483b0559 229 return kFALSE;
230 }
231 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
232 index++ ;
233 if(ierflg != 0){
9859bfc0 234 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ;
483b0559 235 return kFALSE;
236 }
237 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
238 index++ ;
239 if(ierflg != 0){
9859bfc0 240 Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ;
483b0559 241 return kFALSE;
242 }
243 }
244
245 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
246 // depends on it.
247 Double_t p1 = 1.0 ;
248 Double_t p2 = 0.0 ;
249
250 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
251 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
252 gMinuit->SetMaxIterations(5);
253 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
254
255 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
256
257 if(ierflg == 4){ // Minimum not found
9859bfc0 258 Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ;
483b0559 259 return kFALSE ;
260 }
261 for(index = 0; index < nPar; index++){
262 Double_t err ;
263 Double_t val ;
264 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
265 fitparameters[index] = val ;
266 }
267
268 delete toMinuit ;
269 return kTRUE;
270
271}
272
273//____________________________________________________________________________
274void AliEMCALClusterizerv1::GetCalibrationParameters()
275{
88cb7938 276 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
483b0559 277
88cb7938 278 if ( !gime->Digitizer() )
279 gime->LoadDigitizer();
280 AliEMCALDigitizer * dig = gime->Digitizer();
281
ed611565 282 fADCchannelPRE = dig->GetPREchannel() ;
283 fADCpedestalPRE = dig->GetPREpedestal() ;
284
88cb7938 285 fADCchannelECA = dig->GetECAchannel() ;
286 fADCpedestalECA = dig->GetECApedestal();
483b0559 287
88cb7938 288 fADCchannelHCA = dig->GetHCAchannel() ;
289 fADCpedestalHCA = dig->GetHCApedestal();
483b0559 290}
05a92d59 291
483b0559 292//____________________________________________________________________________
293void AliEMCALClusterizerv1::Init()
294{
295 // Make all memory allocations which can not be done in default constructor.
296 // Attach the Clusterizer task to the list of EMCAL tasks
297
88cb7938 298 AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data());
483b0559 299
88cb7938 300 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
05a92d59 301
483b0559 302 fNTowers = geom->GetNZ() * geom->GetNPhi() ;
483b0559 303 if(!gMinuit)
304 gMinuit = new TMinuit(100) ;
305
88cb7938 306 if ( !gime->Clusterizer() )
307 gime->PostClusterizer(this);
483b0559 308}
309
310//____________________________________________________________________________
839828a6 311void AliEMCALClusterizerv1::InitParameters()
312{
88cb7938 313 fNumberOfPREClusters = fNumberOfECAClusters = fNumberOfHCAClusters = 0 ;
ed611565 314 fPREClusteringThreshold = 0.0001; // must be adjusted according to the noise leve set by digitizer
88cb7938 315 fECAClusteringThreshold = 0.0045; // must be adjusted according to the noise leve set by digitizer
316 fHCAClusteringThreshold = 0.001; // must be adjusted according to the noise leve set by digitizer
ed611565 317 fPRELocMaxCut = 0.03 ;
88cb7938 318 fECALocMaxCut = 0.03 ;
319 fHCALocMaxCut = 0.03 ;
839828a6 320
ed611565 321 fPREW0 = 4.0 ;
88cb7938 322 fECAW0 = 4.5 ;
323 fHCAW0 = 4.5 ;
839828a6 324
325 fTimeGate = 1.e-8 ;
326
327 fToUnfold = kFALSE ;
88cb7938 328
329 fRecPointsInRun = 0 ;
330
839828a6 331}
332
333//____________________________________________________________________________
483b0559 334Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const
335{
336 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
337 // = 1 are neighbour
338 // = 2 are not neighbour but do not continue searching
339 // neighbours are defined as digits having at least a common vertex
340 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
341 // which is compared to a digit (d2) not yet in a cluster
342
88cb7938 343 AliEMCALGeometry * geom = AliEMCALGetter::Instance()->EMCALGeometry() ;
483b0559 344
345 Int_t rv = 0 ;
346
347 Int_t relid1[4] ;
348 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
349
350 Int_t relid2[4] ;
351 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
ed611565 352
353 if ( (relid1[0] == relid2[0]) && // inside the same EMCAL Arm
354 (relid1[1]==relid2[1]) ) { // and same tower section
483b0559 355 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
356 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
357
358 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
359 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate))
360 rv = 1 ;
361 }
362 else {
363 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
364 rv = 2; // Difference in row numbers is too large to look further
365 }
366
367 }
368 else {
369
370 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
ed611565 371 rv=0 ;
483b0559 372 }
483b0559 373
ed611565 374 if (gDebug == 2 )
375 Info("AreNeighbours", "neighbours=%d, id1=%d, relid1=%d,%d,%d,%d \n id2=%d, relid2=%d,%d,%d,%d",
376 rv, d1->GetId(), relid1[0], relid1[1], relid1[2], relid1[3],
377 d2->GetId(), relid2[0], relid2[1], relid2[2], relid2[3]) ;
378
483b0559 379 return rv ;
380}
381
88cb7938 382//____________________________________________________________________________
383void AliEMCALClusterizerv1::Unload()
384{
385 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
386 gime->EmcalLoader()->UnloadDigits() ;
387 gime->EmcalLoader()->UnloadRecPoints() ;
388}
483b0559 389
483b0559 390//____________________________________________________________________________
391void AliEMCALClusterizerv1::WriteRecPoints(Int_t event)
392{
393
394 // Creates new branches with given title
395 // fills and writes into TreeR.
396
88cb7938 397 AliEMCALGetter *gime = AliEMCALGetter::Instance() ;
ed611565 398
399 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
88cb7938 400 TObjArray * aECARecPoints = gime->ECARecPoints() ;
401 TObjArray * aHCARecPoints = gime->HCARecPoints() ;
ed611565 402
05a92d59 403 TClonesArray * digits = gime->Digits() ;
88cb7938 404 TTree * treeR = gime->TreeR(); ;
05a92d59 405
05a92d59 406 Int_t index ;
ed611565 407
408 //Evaluate position, dispersion and other RecPoint properties for PRE section
409 for(index = 0; index < aPRERecPoints->GetEntries(); index++)
410 (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->EvalAll(fPREW0,digits) ;
411 aPRERecPoints->Sort() ;
412
413 for(index = 0; index < aPRERecPoints->GetEntries(); index++)
414 (dynamic_cast<AliEMCALRecPoint *>(aPRERecPoints->At(index)))->SetIndexInList(index) ;
415
416 aPRERecPoints->Expand(aPRERecPoints->GetEntriesFast()) ;
05a92d59 417
ed611565 418 //Evaluate position, dispersion and other RecPoint properties for EC section
88cb7938 419 for(index = 0; index < aECARecPoints->GetEntries(); index++)
420 (dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ;
ed611565 421
88cb7938 422 aECARecPoints->Sort() ;
483b0559 423
88cb7938 424 for(index = 0; index < aECARecPoints->GetEntries(); index++)
425 (dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(index)))->SetIndexInList(index) ;
483b0559 426
88cb7938 427 aECARecPoints->Expand(aECARecPoints->GetEntriesFast()) ;
ed611565 428
88cb7938 429 //Evaluate position, dispersion and other RecPoint properties for HCA section
430 for(index = 0; index < aHCARecPoints->GetEntries(); index++)
431 (dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(index)))->EvalAll(fHCAW0,digits) ;
ed611565 432
88cb7938 433 aHCARecPoints->Sort() ;
483b0559 434
88cb7938 435 for(index = 0; index < aHCARecPoints->GetEntries(); index++)
436 (dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(index)))->SetIndexInList(index) ;
483b0559 437
88cb7938 438 aHCARecPoints->Expand(aHCARecPoints->GetEntriesFast()) ;
ed611565 439
483b0559 440 Int_t bufferSize = 32000 ;
9859bfc0 441 Int_t splitlevel = 0 ;
05a92d59 442
ed611565 443 //PRE section branch
444 TBranch * branchPRE = treeR->Branch("EMCALPRERP","TObjArray",&aPRERecPoints,bufferSize,splitlevel);
445 branchPRE->SetTitle(BranchName());
446
447 //EC section branch
88cb7938 448 TBranch * branchECA = treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel);
449 branchECA->SetTitle(BranchName());
ed611565 450
88cb7938 451 //HCA section branch
452 TBranch * branchHCA = treeR->Branch("EMCALHCARP","TObjArray",&aHCARecPoints,bufferSize,splitlevel);
453 branchHCA->SetTitle(BranchName());
839828a6 454
ed611565 455 branchPRE->Fill() ;
88cb7938 456 branchECA->Fill() ;
457 branchHCA->Fill() ;
483b0559 458
88cb7938 459 gime->WriteRecPoints("OVERWRITE");
460 gime->WriteClusterizer("OVERWRITE");
483b0559 461}
462
463//____________________________________________________________________________
464void AliEMCALClusterizerv1::MakeClusters()
465{
466 // Steering method to construct the clusters stored in a list of Reconstructed Points
467 // A cluster is defined as a list of neighbour digits
05a92d59 468
88cb7938 469 AliEMCALGetter * gime = AliEMCALGetter::Instance() ;
ed611565 470
471 AliEMCALGeometry * geom = gime->EMCALGeometry() ;
472
473
88cb7938 474 TObjArray * aPRERecPoints = gime->PRERecPoints() ;
475 TObjArray * aECARecPoints = gime->ECARecPoints() ;
476 TObjArray * aHCARecPoints = gime->HCARecPoints() ;
ed611565 477
478 aPRERecPoints->Delete() ;
88cb7938 479 aECARecPoints->Delete() ;
480 aHCARecPoints->Delete() ;
ed611565 481
05a92d59 482 TClonesArray * digits = gime->Digits() ;
ed611565 483
484 TIter next(digits) ;
485 AliEMCALDigit * digit ;
88cb7938 486 Int_t ndigECA=0, ndigPRE=0, ndigHCA=0 ;
ed611565 487
88cb7938 488 // count the number of digits in ECA section
ed611565 489 while ( (digit = dynamic_cast<AliEMCALDigit *>(next())) ) { // scan over the list of digits
88cb7938 490 if (geom->IsInECA(digit->GetId()))
491 ndigECA++ ;
ed611565 492 else if (geom->IsInPRE(digit->GetId()))
493 ndigPRE++;
88cb7938 494 else if (geom->IsInHCA(digit->GetId()))
495 ndigHCA++;
ed611565 496 else {
497 Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ;
498 abort() ;
499 }
500 }
501
88cb7938 502 // add amplitude of PRE and ECA sections
503 Int_t digECA ;
504 for (digECA = 0 ; digECA < ndigECA ; digECA++) {
505 digit = dynamic_cast<AliEMCALDigit *>(digits->At(digECA)) ;
ed611565 506 Int_t digPRE ;
88cb7938 507 for (digPRE = ndigECA ; digPRE < ndigECA+ndigPRE ; digPRE++) {
ed611565 508 AliEMCALDigit * digitPRE = dynamic_cast<AliEMCALDigit *>(digits->At(digPRE)) ;
509 if ( geom->AreInSameTower(digit->GetId(), digitPRE->GetId()) ){
510 Float_t amp = static_cast<Float_t>(digit->GetAmp()) + geom->GetSummationFraction() * static_cast<Float_t>(digitPRE->GetAmp()) + 0.5 ;
511 digit->SetAmp(static_cast<Int_t>(amp)) ;
512 break ;
513 }
514 }
515 if (gDebug)
516 Info("MakeClusters","id = %d amp = %d", digit->GetId(), digit->GetAmp()) ;
517 }
518
483b0559 519 TClonesArray * digitsC = dynamic_cast<TClonesArray*>(digits->Clone()) ;
520
521
522 // Clusterization starts
523
524 TIter nextdigit(digitsC) ;
88cb7938 525 Bool_t notremovedECA = kTRUE, notremovedPRE = kTRUE ;
483b0559 526
527 while ( (digit = dynamic_cast<AliEMCALDigit *>(nextdigit())) ) { // scan over the list of digitsC
528 AliEMCALRecPoint * clu = 0 ;
529
88cb7938 530 TArrayI clusterPREdigitslist(50), clusterECAdigitslist(50), clusterHCAdigitslist(50);
ed611565 531
88cb7938 532 Bool_t inPRE = kFALSE, inECA = kFALSE, inHCA = kFALSE ;
ed611565 533 if( geom->IsInPRE(digit->GetId()) ) {
534 inPRE = kTRUE ;
535 }
88cb7938 536 else if( geom->IsInECA(digit->GetId()) ) {
537 inECA = kTRUE ;
ed611565 538 }
88cb7938 539 else if( geom->IsInHCA(digit->GetId()) ) {
540 inHCA = kTRUE ;
ed611565 541 }
542
543 if (gDebug == 2) {
544 if (inPRE)
545 Info("MakeClusters","id = %d, ene = %f , thre = %f ",
546 digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
88cb7938 547 if (inECA)
ed611565 548 Info("MakeClusters","id = %d, ene = %f , thre = %f",
88cb7938 549 digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ;
550 if (inHCA)
ed611565 551 Info("MakeClusters","id = %d, ene = %f , thre = %f",
88cb7938 552 digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold ) ;
ed611565 553 }
554
555 if ( (inPRE && (Calibrate(digit->GetAmp(), 0) > fPREClusteringThreshold )) ||
88cb7938 556 (inECA && (Calibrate(digit->GetAmp(), 1) > fECAClusteringThreshold )) ||
557 (inHCA && (Calibrate(digit->GetAmp(), 2) > fHCAClusteringThreshold )) ) {
483b0559 558
88cb7938 559 Int_t iDigitInPRECluster = 0, iDigitInECACluster = 0, iDigitInHCACluster = 0;
ed611565 560 Int_t where ; // PRE = 0, ECAl = 1, HCAL = 2
561
562 // Find the seed in each of the section ECAL/PRE/HCAL
563
88cb7938 564 if( geom->IsInECA(digit->GetId()) ) {
ed611565 565 where = 1 ; // to tell we are in ECAL
483b0559 566 // start a new Tower RecPoint
88cb7938 567 if(fNumberOfECAClusters >= aECARecPoints->GetSize())
568 aECARecPoints->Expand(2*fNumberOfECAClusters+1) ;
37b680d1 569 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
88cb7938 570 rp->SetECA() ;
571 aECARecPoints->AddAt(rp, fNumberOfECAClusters) ;
572 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aECARecPoints->At(fNumberOfECAClusters)) ;
573 fNumberOfECAClusters++ ;
ed611565 574 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)) ;
88cb7938 575 clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ;
576 iDigitInECACluster++ ;
483b0559 577 digitsC->Remove(digit) ;
ed611565 578 if (gDebug == 2 )
88cb7938 579 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 1), fECAClusteringThreshold) ;
483b0559 580
ed611565 581 }
582 else if( geom->IsInPRE(digit->GetId()) ) {
583 where = 0 ; // to tell we are in PRE
483b0559 584 // start a new Pre Shower cluster
ed611565 585 if(fNumberOfPREClusters >= aPRERecPoints->GetSize())
586 aPRERecPoints->Expand(2*fNumberOfPREClusters+1);
587 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
588 rp->SetPRE() ;
589 aPRERecPoints->AddAt(rp, fNumberOfPREClusters) ;
590 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(fNumberOfPREClusters)) ;
591 fNumberOfPREClusters++ ;
592 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
593 clusterPREdigitslist[iDigitInPRECluster] = digit->GetIndexInList() ;
594 iDigitInPRECluster++ ;
595 digitsC->Remove(digit) ;
596 if (gDebug == 2 )
597 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ;
483b0559 598
483b0559 599 nextdigit.Reset() ;
ed611565 600
88cb7938 601 // Here we remove remaining ECA digits, which cannot make a cluster
483b0559 602
88cb7938 603 if( notremovedECA ) {
483b0559 604 while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
88cb7938 605 if( geom->IsInECA(digit->GetId()) )
483b0559 606 digitsC->Remove(digit) ;
607 else
608 break ;
609 }
88cb7938 610 notremovedECA = kFALSE ;
483b0559 611 }
ed611565 612
613 }
88cb7938 614 else if( geom->IsInHCA(digit->GetId()) ) {
ed611565 615 where = 2 ; // to tell we are in HCAL
616 // start a new HCAL cluster
88cb7938 617 if(fNumberOfHCAClusters >= aHCARecPoints->GetSize())
618 aHCARecPoints->Expand(2*fNumberOfHCAClusters+1);
ed611565 619 AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ;
88cb7938 620 rp->SetHCA() ;
621 aHCARecPoints->AddAt(rp, fNumberOfHCAClusters) ;
622 clu = dynamic_cast<AliEMCALTowerRecPoint *>(aHCARecPoints->At(fNumberOfHCAClusters)) ;
623 fNumberOfHCAClusters++ ;
ed611565 624 clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where));
88cb7938 625 clusterHCAdigitslist[iDigitInHCACluster] = digit->GetIndexInList() ;
626 iDigitInHCACluster++ ;
ed611565 627 digitsC->Remove(digit) ;
628 if (gDebug == 2 )
88cb7938 629 Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCAClusteringThreshold) ;
ed611565 630
631 nextdigit.Reset() ;
632
633 // Here we remove remaining PRE digits, which cannot make a cluster
483b0559 634
ed611565 635 if( notremovedPRE ) {
636 while( ( digit = dynamic_cast<AliEMCALDigit *>(nextdigit()) ) ) {
637 if( geom->IsInPRE(digit->GetId()) )
638 digitsC->Remove(digit) ;
639 else
640 break ;
641 }
642 notremovedPRE = kFALSE ;
643 }
644 }
483b0559 645
646 nextdigit.Reset() ;
647
648 AliEMCALDigit * digitN ;
ed611565 649 Int_t index = 0 ;
650
651 // Do the Clustering in each of the three section ECAL/PRE/HCAL
652
88cb7938 653 while (index < iDigitInECACluster){ // scan over digits already in cluster
654 digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]) ;
483b0559 655 index++ ;
656 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
657 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
ed611565 658 // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
483b0559 659 switch (ineb ) {
660 case 0 : // not a neighbour
661 break ;
662 case 1 : // are neighbours
ed611565 663 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 1) ) ;
88cb7938 664 clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ;
665 iDigitInECACluster++ ;
483b0559 666 digitsC->Remove(digitN) ;
667 break ;
668 case 2 : // too far from each other
ed611565 669 goto endofloop1;
483b0559 670 } // switch
671
672 } // while digitN
673
ed611565 674 endofloop1: ;
483b0559 675 nextdigit.Reset() ;
88cb7938 676 } // loop over ECA cluster
ed611565 677
678 index = 0 ;
679 while (index < iDigitInPRECluster){ // scan over digits already in cluster
680 digit = (AliEMCALDigit*)digits->At(clusterPREdigitslist[index]) ;
681 index++ ;
682 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
683 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
684 // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
685 switch (ineb ) {
686 case 0 : // not a neighbour
687 break ;
688 case 1 : // are neighbours
689 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 0) ) ;
690 clusterPREdigitslist[iDigitInPRECluster] = digitN->GetIndexInList() ;
691 iDigitInPRECluster++ ;
692 digitsC->Remove(digitN) ;
693 break ;
694 case 2 : // too far from each other
695 goto endofloop2;
696 } // switch
697
698 } // while digitN
699
700 endofloop2: ;
701 nextdigit.Reset() ;
702 } // loop over PRE cluster
703
704 index = 0 ;
88cb7938 705 while (index < iDigitInHCACluster){ // scan over digits already in cluster
706 digit = (AliEMCALDigit*)digits->At(clusterHCAdigitslist[index]) ;
ed611565 707 index++ ;
708 while ( (digitN = (AliEMCALDigit *)nextdigit()) ) { // scan over the reduced list of digits
709 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
710 //Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ;
711 switch (ineb ) {
712 case 0 : // not a neighbour
713 break ;
714 case 1 : // are neighbours
715 clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 2) ) ;
88cb7938 716 clusterHCAdigitslist[iDigitInHCACluster] = digitN->GetIndexInList() ;
717 iDigitInHCACluster++ ;
ed611565 718 digitsC->Remove(digitN) ;
719 break ;
720 case 2 : // too far from each other
721 goto endofloop3;
722 } // switch
723 } // while digitN
724
725 endofloop3: ;
726 nextdigit.Reset() ;
88cb7938 727 } // loop over HCA cluster
ed611565 728
9859bfc0 729 } // energy theshold
730 } // while digit
483b0559 731 delete digitsC ;
483b0559 732}
733
734//____________________________________________________________________________
735void AliEMCALClusterizerv1::MakeUnfolding()
736{
737 Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ;
9859bfc0 738
483b0559 739}
740
741//____________________________________________________________________________
742Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r)
743{
744 // Shape of the shower (see EMCAL TDR)
745 // If you change this function, change also the gradient evaluation in ChiSquare()
746
747 Double_t r4 = r*r*r*r ;
748 Double_t r295 = TMath::Power(r, 2.95) ;
749 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
750 return shape ;
751}
752
753//____________________________________________________________________________
754void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower,
755 Int_t nMax,
a0636361 756 AliEMCALDigit ** maxAt,
483b0559 757 Float_t * maxAtEnergy)
758{
759 // Performs the unfolding of a cluster with nMax overlapping showers
483b0559 760
9859bfc0 761 Fatal("UnfoldCluster", "--> Unfolding not implemented") ;
173558f2 762
483b0559 763}
764
765//_____________________________________________________________________________
766void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
767{
768 // Calculates the Chi square for the cluster unfolding minimization
769 // Number of parameters, Gradient, Chi squared, parameters, what to do
483b0559 770
9859bfc0 771 ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ;
483b0559 772}
483b0559 773//____________________________________________________________________________
774void AliEMCALClusterizerv1::Print(Option_t * option)const
775{
776 // Print clusterizer parameters
777
9859bfc0 778 TString message("\n") ;
779
483b0559 780 if( strcmp(GetName(), "") !=0 ){
781
782 // Print parameters
783
784 TString taskName(GetName()) ;
785 taskName.ReplaceAll(Version(), "") ;
9859bfc0 786
787 message += "--------------- " ;
788 message += taskName.Data() ;
789 message += " " ;
790 message += GetTitle() ;
791 message += "-----------\n" ;
792 message += "Clusterizing digits from the file: " ;
793 message += taskName.Data() ;
794 message += "\n Branch: " ;
795 message += GetName() ;
f6eaf97a 796 message += "\n Pre Shower Clustering threshold = " ;
ed611565 797 message += fPREClusteringThreshold ;
f6eaf97a 798 message += "\n Pre Shower Local Maximum cut = " ;
ed611565 799 message += fPRELocMaxCut ;
f6eaf97a 800 message += "\n Pre Shower Logarothmic weight = " ;
ed611565 801 message += fPREW0 ;
88cb7938 802 message += "\n ECA Clustering threshold = " ;
803 message += fECAClusteringThreshold ;
804 message += "\n ECA Local Maximum cut = " ;
805 message += fECALocMaxCut ;
806 message += "\n ECA Logarothmic weight = " ;
807 message += fECAW0 ;
ed611565 808 message += "\n Pre Shower Clustering threshold = " ;
88cb7938 809 message += fHCAClusteringThreshold ;
810 message += "\n HCA Local Maximum cut = " ;
811 message += fHCALocMaxCut ;
812 message += "\n HCA Logarothmic weight = " ;
813 message += fHCAW0 ;
483b0559 814 if(fToUnfold)
9859bfc0 815 message +="\nUnfolding on\n" ;
483b0559 816 else
9859bfc0 817 message += "\nUnfolding off\n";
483b0559 818
9859bfc0 819 message += "------------------------------------------------------------------" ;
483b0559 820 }
821 else
9859bfc0 822 message += "AliEMCALClusterizerv1 not initialized " ;
823
824 Info("Print", message.Data() ) ;
483b0559 825}
173558f2 826
483b0559 827//____________________________________________________________________________
828void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option)
829{
830 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
831
88cb7938 832 TObjArray * aPRERecPoints = AliEMCALGetter::Instance()->PRERecPoints() ;
833 TObjArray * aECARecPoints = AliEMCALGetter::Instance()->ECARecPoints() ;
834 TObjArray * aHCARecPoints = AliEMCALGetter::Instance()->HCARecPoints() ;
173558f2 835
ed611565 836 Info("PrintRecPoints", "Clusterization result:") ;
837
838 printf("event # %d\n", gAlice->GetEvNumber() ) ;
88cb7938 839 printf(" Found %d PRE SHOWER RecPoints, %d ECA Rec Points and %d HCA Rec Points\n ",
840 aPRERecPoints->GetEntriesFast(), aECARecPoints->GetEntriesFast(), aHCARecPoints->GetEntriesFast() ) ;
483b0559 841
ed611565 842 fRecPointsInRun += aPRERecPoints->GetEntriesFast() ;
88cb7938 843 fRecPointsInRun += aECARecPoints->GetEntriesFast() ;
844 fRecPointsInRun += aHCARecPoints->GetEntriesFast() ;
11f9c5ff 845
483b0559 846 if(strstr(option,"all")) {
847
ed611565 848 //Pre shower recPoints
849
850 printf("-----------------------------------------------------------------------\n") ;
851 printf("Clusters in PRE section\n") ;
852 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
853
483b0559 854 Int_t index ;
ed611565 855
856 for (index = 0 ; index < aPRERecPoints->GetEntries() ; index++) {
857 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint *>(aPRERecPoints->At(index)) ;
858 TVector3 globalpos;
859 rp->GetGlobalPosition(globalpos);
860 TVector3 localpos;
861 rp->GetLocalPosition(localpos);
862 Float_t lambda[2];
863 rp->GetElipsAxis(lambda);
864 Int_t * primaries;
865 Int_t nprimaries;
866 primaries = rp->GetPrimaries(nprimaries);
867 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
868 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
869 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
870 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
871 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
872 printf("%d ", primaries[iprimary] ) ;
873 }
874 }
875
876 printf("\n-----------------------------------------------------------------------\n") ;
877 printf("Clusters in ECAL section\n") ;
878 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
879
88cb7938 880 for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) {
881 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aECARecPoints->At(index)) ;
483b0559 882 TVector3 globalpos;
883 rp->GetGlobalPosition(globalpos);
ed611565 884 TVector3 localpos;
885 rp->GetLocalPosition(localpos);
483b0559 886 Float_t lambda[2];
887 rp->GetElipsAxis(lambda);
888 Int_t * primaries;
889 Int_t nprimaries;
890 primaries = rp->GetPrimaries(nprimaries);
ed611565 891 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
892 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
893 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
894 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
9859bfc0 895 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
ed611565 896 printf("%d ", primaries[iprimary] ) ;
9859bfc0 897 }
483b0559 898 }
899
ed611565 900 printf("\n-----------------------------------------------------------------------\n") ;
901 printf("Clusters in HCAL section\n") ;
902 printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
483b0559 903
88cb7938 904 for (index = 0 ; index < aHCARecPoints->GetEntries() ; index++) {
905 AliEMCALTowerRecPoint * rp = dynamic_cast<AliEMCALTowerRecPoint * >(aHCARecPoints->At(index)) ;
483b0559 906 TVector3 globalpos;
907 rp->GetGlobalPosition(globalpos);
ed611565 908 TVector3 localpos;
909 rp->GetLocalPosition(localpos);
483b0559 910 Float_t lambda[2];
911 rp->GetElipsAxis(lambda);
ed611565 912 Int_t * primaries;
483b0559 913 Int_t nprimaries;
914 primaries = rp->GetPrimaries(nprimaries);
ed611565 915 printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
916 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(),
917 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
918 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
9859bfc0 919 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
ed611565 920 printf("%d ", primaries[iprimary] ) ;
921 }
483b0559 922 }
923
ed611565 924 printf("\n-----------------------------------------------------------------------\n");
483b0559 925 }
926}