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