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