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