High/Low and time calibrations added
[u/mrichter/AliRoot.git] / PHOS / AliPHOSClusterizerv1.cxx
CommitLineData
d15a28e7 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 **************************************************************************/
15
b2a60966 16/* $Id$ */
17
702ab87e 18/* History of cvs commits:
19 *
20 * $Log$
a42611a9 21 * Revision 1.117 2007/10/18 08:42:05 kharlov
22 * Bad channels cleaned before clusterization
23 *
454c62a4 24 * Revision 1.116 2007/10/01 20:24:08 kharlov
25 * Memory leaks fixed
26 *
e68222ce 27 * Revision 1.115 2007/09/26 14:22:17 cvetan
28 * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods.
29 *
d76c31f4 30 * Revision 1.114 2007/09/06 16:06:44 kharlov
31 * Absence of sorting results in loose of all unfolded clusters
32 *
6d0a5c2e 33 * Revision 1.113 2007/08/28 12:55:07 policheh
34 * Loaders removed from the reconstruction code (C.Cheshkov)
35 *
9a2cdbdf 36 * Revision 1.112 2007/08/22 09:20:50 hristov
37 * Updated QA classes (Yves)
38 *
a5fa6165 39 * Revision 1.111 2007/08/08 12:11:28 kharlov
40 * Protection against uninitialized fQADM
41 *
7ae07b77 42 * Revision 1.110 2007/08/07 14:16:00 kharlov
43 * Quality assurance added (Yves Schutz)
44 *
ddd1a39c 45 * Revision 1.109 2007/07/24 17:20:35 policheh
46 * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
47 * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
48 *
3799bcb5 49 * Revision 1.108 2007/06/18 07:00:51 kharlov
50 * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
51 *
cd79ec76 52 * Revision 1.107 2007/05/25 14:12:26 policheh
53 * Local to tracking CS transformation added for CPV rec. points
54 *
76e42795 55 * Revision 1.106 2007/05/24 13:01:22 policheh
56 * Local to tracking CS transformation invoked for each EMC rec.point
57 *
788ab371 58 * Revision 1.105 2007/05/02 13:41:22 kharlov
59 * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
60 *
f36dfa9a 61 * Revision 1.104 2007/04/27 16:55:53 kharlov
62 * Calibration stops if PHOS CDB objects do not exist
63 *
76d78b7a 64 * Revision 1.103 2007/04/11 11:55:45 policheh
65 * SetDistancesToBadChannels() added.
66 *
af0b0570 67 * Revision 1.102 2007/03/28 19:18:15 kharlov
68 * RecPoints recalculation in TSM removed
69 *
8c1fb709 70 * Revision 1.101 2007/03/06 06:51:27 kharlov
71 * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
72 *
a978bcea 73 * Revision 1.100 2007/01/10 11:07:26 kharlov
74 * Raw digits writing to file (B.Polichtchouk)
75 *
cbddd97f 76 * Revision 1.99 2006/11/07 16:49:51 kharlov
77 * Corrections for next event switch in case of raw data (B.Polichtchouk)
78 *
e1c29d81 79 * Revision 1.98 2006/10/27 17:14:27 kharlov
80 * Introduce AliDebug and AliLog (B.Polichtchouk)
81 *
e1afea68 82 * Revision 1.97 2006/08/29 11:41:19 kharlov
83 * Missing implementation of ctors and = operator are added
84 *
f1133801 85 * Revision 1.96 2006/08/25 16:56:30 kharlov
86 * Compliance with Effective C++
87 *
0378398c 88 * Revision 1.95 2006/08/11 12:36:26 cvetan
89 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
90 *
52c5f046 91 * Revision 1.94 2006/08/07 12:27:49 hristov
92 * Removing obsolete code which affected the event numbering scheme
93 *
4d4750bd 94 * Revision 1.93 2006/08/01 12:20:17 cvetan
95 * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
96 *
f5eaa851 97 * Revision 1.92 2006/04/29 20:26:46 hristov
98 * Separate EMC and CPV calibration (Yu.Kharlov)
99 *
e95226ae 100 * Revision 1.91 2006/04/22 10:30:17 hristov
101 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
102 *
27a73a5d 103 * Revision 1.90 2006/04/11 15:22:59 hristov
104 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
105 *
28871337 106 * Revision 1.89 2006/03/13 14:05:42 kharlov
107 * Calibration objects for EMC and CPV
108 *
fc6706cb 109 * Revision 1.88 2006/01/11 08:54:52 hristov
110 * Additional protection in case no calibration entry was found
111 *
e521c5a0 112 * Revision 1.87 2005/11/22 08:46:43 kharlov
113 * Updated to new CDB (Boris Polichtchouk)
114 *
ef629623 115 * Revision 1.86 2005/11/14 21:52:43 hristov
116 * Coding conventions
117 *
fdf65bb5 118 * Revision 1.85 2005/09/27 16:08:08 hristov
119 * New version of CDB storage framework (A.Colla)
120 *
9e1ceb13 121 * Revision 1.84 2005/09/21 10:02:47 kharlov
122 * Reading calibration from CDB (Boris Polichtchouk)
123 *
8d6d4fb5 124 * Revision 1.82 2005/09/02 15:43:13 kharlov
125 * Add comments in GetCalibrationParameters and Calibrate
126 *
5c149aeb 127 * Revision 1.81 2005/09/02 14:32:07 kharlov
128 * Calibration of raw data
129 *
44ae287e 130 * Revision 1.80 2005/08/24 15:31:36 kharlov
131 * Setting raw digits flag
132 *
9735c615 133 * Revision 1.79 2005/07/25 15:53:53 kharlov
134 * Read raw data
135 *
3cf4f75f 136 * Revision 1.78 2005/05/28 14:19:04 schutz
137 * Compilation warnings fixed by T.P.
138 *
702ab87e 139 */
140
9a1398dd 141//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
d15a28e7 142//////////////////////////////////////////////////////////////////////////////
9a1398dd 143// Clusterization class. Performs clusterization (collects neighbouring active cells) and
a4e98857 144// unfolds the clusters having several local maxima.
145// Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
9a1398dd 146// PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
f035f6ce 147// parameters including input digits branch title, thresholds etc.)
a4e98857 148// This TTask is normally called from Reconstructioner, but can as well be used in
149// standalone mode.
150// Use Case:
9a2cdbdf 151// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
152// root [1] cl->Digits2Clusters(digitsTree,clusterTree)
153// //finds RecPoints in the current event
a4e98857 154// root [2] cl->SetDigitsBranch("digits2")
f035f6ce 155// //sets another title for Digitis (input) branch
a4e98857 156// root [3] cl->SetRecPointsBranch("recp2")
f035f6ce 157// //sets another title four output branches
a4e98857 158// root [4] cl->SetEmcLocalMaxCut(0.03)
9a1398dd 159// //set clusterization parameters
ed4205d8 160
d15a28e7 161// --- ROOT system ---
162
163#include "TMath.h"
9a1398dd 164#include "TMinuit.h"
165#include "TTree.h"
9a1398dd 166#include "TBenchmark.h"
9a2cdbdf 167#include "TClonesArray.h"
d15a28e7 168
169// --- Standard library ---
170
d15a28e7 171// --- AliRoot header files ---
351dd634 172#include "AliLog.h"
9a2cdbdf 173#include "AliConfig.h"
e957fea8 174#include "AliPHOSGeometry.h"
d15a28e7 175#include "AliPHOSClusterizerv1.h"
88cb7938 176#include "AliPHOSEmcRecPoint.h"
9a1398dd 177#include "AliPHOSCpvRecPoint.h"
d15a28e7 178#include "AliPHOSDigit.h"
9a1398dd 179#include "AliPHOSDigitizer.h"
ba898748 180#include "AliPHOSCalibrationDB.h"
9e1ceb13 181#include "AliCDBManager.h"
eaa5012b 182#include "AliCDBStorage.h"
fdf65bb5 183#include "AliCDBEntry.h"
3799bcb5 184#include "AliPHOSRecoParam.h"
9a2cdbdf 185#include "AliPHOSCalibData.h"
186#include "AliPHOSReconstructor.h"
d15a28e7 187
188ClassImp(AliPHOSClusterizerv1)
f035f6ce 189
d15a28e7 190//____________________________________________________________________________
0378398c 191AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
192 AliPHOSClusterizer(),
193 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
194 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
e68222ce 195 fADCchanelEmc(0), fADCpedestalEmc(0),
0378398c 196 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
197 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
198 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
9a2cdbdf 199 fW0CPV(0), fEmcTimeGate(0),
200 fIsOldRCUFormat(0)
d15a28e7 201{
f035f6ce 202 // default ctor (to be used mainly by Streamer)
f035f6ce 203
8d0f3f77 204 InitParameters() ;
92f521a9 205 fDefaultInit = kTRUE ;
9a1398dd 206}
7b7c1533 207
9a1398dd 208//____________________________________________________________________________
9a2cdbdf 209AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
210 AliPHOSClusterizer(geom),
0378398c 211 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
212 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
e68222ce 213 fADCchanelEmc(0), fADCpedestalEmc(0),
0378398c 214 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
215 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
216 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
9a2cdbdf 217 fW0CPV(0), fEmcTimeGate(0),
218 fIsOldRCUFormat(0)
9a1398dd 219{
a4e98857 220 // ctor with the indication of the file where header Tree and digits Tree are stored
9a1398dd 221
8d0f3f77 222 InitParameters() ;
2bd5457f 223 Init() ;
92f521a9 224 fDefaultInit = kFALSE ;
9a1398dd 225}
fc12304f 226
9688c1dd 227//____________________________________________________________________________
228 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
229{
0bc3b8ed 230 // dtor
231
9688c1dd 232}
9a2cdbdf 233
9a1398dd 234//____________________________________________________________________________
e95226ae 235Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
88cb7938 236{
e95226ae 237 // Convert EMC measured amplitude into real energy.
5c149aeb 238 // Calibration parameters are taken from calibration data base for raw data,
239 // or from digitizer parameters for simulated data.
240
e68222ce 241 if(fgCalibData){
44ae287e 242 Int_t relId[4];
9a2cdbdf 243 fGeom->AbsToRelNumbering(absId,relId) ;
44ae287e 244 Int_t module = relId[0];
8d6d4fb5 245 Int_t column = relId[3];
eaa5012b 246 Int_t row = relId[2];
e95226ae 247 if(absId <= fEmcCrystals) { // this is EMC
e68222ce 248 fADCchanelEmc = fgCalibData->GetADCchannelEmc (module,column,row);
e95226ae 249 return amp*fADCchanelEmc ;
44ae287e 250 }
e95226ae 251 }
252 else{ //simulation
253 if(absId <= fEmcCrystals) // this is EMC
254 return fADCpedestalEmc + amp*fADCchanelEmc ;
255 }
256 return 0;
257}
258
259//____________________________________________________________________________
260Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
261{
262 // Convert digitized CPV amplitude into charge.
263 // Calibration parameters are taken from calibration data base for raw data,
264 // or from digitizer parameters for simulated data.
265
e68222ce 266 if(fgCalibData){
e95226ae 267 Int_t relId[4];
9a2cdbdf 268 fGeom->AbsToRelNumbering(absId,relId) ;
e95226ae 269 Int_t module = relId[0];
270 Int_t column = relId[3];
271 Int_t row = relId[2];
272 if(absId > fEmcCrystals) { // this is CPV
e68222ce 273 fADCchanelCpv = fgCalibData->GetADCchannelCpv (module,column,row);
274 fADCpedestalCpv = fgCalibData->GetADCpedestalCpv(module,column,row);
fc6706cb 275 return fADCpedestalCpv + amp*fADCchanelCpv ;
276 }
ba898748 277 }
278 else{ //simulation
e95226ae 279 if(absId > fEmcCrystals) // this is CPV
ba898748 280 return fADCpedestalCpv+ amp*fADCchanelCpv ;
281 }
e95226ae 282 return 0;
3758d9fc 283}
548f0134 284
3758d9fc 285//____________________________________________________________________________
9a2cdbdf 286void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
a4e98857 287{
9a2cdbdf 288 // Steering method to perform clusterization for one event
289 // The input is the tree with digits.
290 // The output is the tree with clusters.
9a1398dd 291
292 if(strstr(option,"tim"))
7d493c2c 293 gBenchmark->Start("PHOSClusterizer");
9a1398dd 294
eabde521 295 if(strstr(option,"print")) {
88cb7938 296 Print() ;
eabde521 297 return ;
298 }
bed9e3fb 299
9a2cdbdf 300 MakeClusters() ;
9a1398dd 301
9a2cdbdf 302 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
303 fEMCRecPoints->GetEntries()));
304 if(AliLog::GetGlobalDebugLevel()>1)
305 fEMCRecPoints->Print();
88cb7938 306
9a2cdbdf 307 if(fToUnfold)
308 MakeUnfolding();
ddd1a39c 309
9a2cdbdf 310 WriteRecPoints();
7b7c1533 311
9a2cdbdf 312 if(strstr(option,"deb"))
313 PrintRecPoints(option) ;
94de8339 314
9a1398dd 315 if(strstr(option,"tim")){
316 gBenchmark->Stop("PHOSClusterizer");
9a2cdbdf 317 AliInfo(Form("took %f seconds for Clusterizing\n",
318 gBenchmark->GetCpuTime("PHOSClusterizer")));
88cb7938 319 }
9a1398dd 320}
321
322//____________________________________________________________________________
a0636361 323Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
d1de15f5 324 Int_t nPar, Float_t * fitparameters) const
9a1398dd 325{
326 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
a4e98857 327 // The initial values for fitting procedure are set equal to the positions of local maxima.
f035f6ce 328 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
9a1398dd 329
88cb7938 330
9a1398dd 331 gMinuit->mncler(); // Reset Minuit's list of paramters
332 gMinuit->SetPrintLevel(-1) ; // No Printout
333 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
334 // To set the address of the minimization function
335
336 TList * toMinuit = new TList();
337 toMinuit->AddAt(emcRP,0) ;
9a2cdbdf 338 toMinuit->AddAt(fDigitsArr,1) ;
339 toMinuit->AddAt(fGeom,2) ;
9a1398dd 340
341 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
342
343 // filling initial values for fit parameters
344 AliPHOSDigit * digit ;
345
346 Int_t ierflg = 0;
347 Int_t index = 0 ;
348 Int_t nDigits = (Int_t) nPar / 3 ;
349
350 Int_t iDigit ;
351
9a1398dd 352 for(iDigit = 0; iDigit < nDigits; iDigit++){
a0636361 353 digit = maxAt[iDigit];
9a1398dd 354
355 Int_t relid[4] ;
7b7c1533 356 Float_t x = 0.;
357 Float_t z = 0.;
9a2cdbdf 358 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
359 fGeom->RelPosInModule(relid, x, z) ;
9a1398dd 360
361 Float_t energy = maxAtEnergy[iDigit] ;
362
363 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
364 index++ ;
365 if(ierflg != 0){
21cd0c07 366 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
9a1398dd 367 return kFALSE;
368 }
369 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
370 index++ ;
371 if(ierflg != 0){
21cd0c07 372 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
9a1398dd 373 return kFALSE;
374 }
375 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
376 index++ ;
377 if(ierflg != 0){
21cd0c07 378 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
9a1398dd 379 return kFALSE;
380 }
381 }
382
383 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
384 // depends on it.
385 Double_t p1 = 1.0 ;
386 Double_t p2 = 0.0 ;
387
388 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
389 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
390 gMinuit->SetMaxIterations(5);
391 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
392
393 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
394
395 if(ierflg == 4){ // Minimum not found
21cd0c07 396 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
9a1398dd 397 return kFALSE ;
398 }
399 for(index = 0; index < nPar; index++){
400 Double_t err ;
401 Double_t val ;
402 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
403 fitparameters[index] = val ;
404 }
405
406 delete toMinuit ;
407 return kTRUE;
c3c187e5 408
d15a28e7 409}
410
3758d9fc 411//____________________________________________________________________________
412void AliPHOSClusterizerv1::GetCalibrationParameters()
413{
5c149aeb 414 // Set calibration parameters:
eaa5012b 415 // if calibration database exists, they are read from database,
76d78b7a 416 // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
5c149aeb 417 //
ef629623 418 // It is a user responsilibity to open CDB before reconstruction, for example:
419 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
fbf811ec 420
e68222ce 421 if (!fgCalibData)
422 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
423 if (fgCalibData->GetCalibDataEmc() == 0)
f36dfa9a 424 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
e68222ce 425 if (fgCalibData->GetCalibDataCpv() == 0)
f36dfa9a 426 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
427
3758d9fc 428}
548f0134 429
d15a28e7 430//____________________________________________________________________________
a4e98857 431void AliPHOSClusterizerv1::Init()
432{
2bd5457f 433 // Make all memory allocations which can not be done in default constructor.
434 // Attach the Clusterizer task to the list of PHOS tasks
ba898748 435
9a2cdbdf 436 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
3758d9fc 437
7b7c1533 438 if(!gMinuit)
88cb7938 439 gMinuit = new TMinuit(100);
fbf811ec 440
e68222ce 441 GetCalibrationParameters() ;
442
9a1398dd 443}
7b7c1533 444
8d0f3f77 445//____________________________________________________________________________
446void AliPHOSClusterizerv1::InitParameters()
447{
448
449 fNumberOfCpvClusters = 0 ;
450 fNumberOfEmcClusters = 0 ;
3799bcb5 451
452 const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
453 if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
454
455 const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv();
456 if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
457
458 fCpvClusteringThreshold = parCpv->GetClusteringThreshold();
459 fEmcClusteringThreshold = parEmc->GetClusteringThreshold();
8d0f3f77 460
3799bcb5 461 fEmcLocMaxCut = parEmc->GetLocalMaxCut();
462 fCpvLocMaxCut = parCpv->GetLocalMaxCut();
ba898748 463
3799bcb5 464 fEmcMinE = parEmc->GetMinE();
465 fCpvMinE = parCpv->GetMinE();
ba898748 466
3799bcb5 467 fW0 = parEmc->GetLogWeight();
468 fW0CPV = parCpv->GetLogWeight();
8d0f3f77 469
454c62a4 470 fEmcTimeGate = 1.e-6 ;
8d0f3f77 471
a42611a9 472 fToUnfold = parEmc->ToUnfold() ;
88cb7938 473
ba898748 474 fWrite = kTRUE ;
475
f5eaa851 476 fIsOldRCUFormat = kFALSE;
8d0f3f77 477}
478
9a1398dd 479//____________________________________________________________________________
480Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
d15a28e7 481{
b2a60966 482 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
483 // = 1 are neighbour
484 // = 2 are not neighbour but do not continue searching
a4e98857 485 // neighbours are defined as digits having at least a common vertex
b2a60966 486 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
487 // which is compared to a digit (d2) not yet in a cluster
488
d15a28e7 489 Int_t rv = 0 ;
490
d15a28e7 491 Int_t relid1[4] ;
9a2cdbdf 492 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
d15a28e7 493
494 Int_t relid2[4] ;
9a2cdbdf 495 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
d15a28e7 496
9688c1dd 497 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
92862013 498 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
499 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
d15a28e7 500
454c62a4 501 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
502 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
9688c1dd 503 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
d15a28e7 504 rv = 1 ;
505 }
506 else {
507 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
88cb7938 508 rv = 2; // Difference in row numbers is too large to look further
d15a28e7 509 }
510
511 }
512 else {
513
9688c1dd 514 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
d15a28e7 515 rv=2 ;
516
517 }
d72dfbc3 518
d15a28e7 519 return rv ;
520}
ba898748 521//____________________________________________________________________________
e95226ae 522void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
523{
454c62a4 524 // Remove digits with amplitudes below threshold.
525 // remove digits in bad channels
e95226ae 526
454c62a4 527 Bool_t isBadMap = 0 ;
528 if(fgCalibData->GetNumOfEmcBadChannels()){
529 isBadMap=1 ;
530 }
531
532 Int_t inBadList=0 ;
ba898748 533 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
534 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
e95226ae 535 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
454c62a4 536 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) ){
ba898748 537 digits->RemoveAt(i) ;
454c62a4 538 continue ;
539 }
540 if(isBadMap){ //check bad map now
541 Int_t relid[4] ;
542 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
543 if(fgCalibData->IsBadChannelEmc(relid[0],relid[2],relid[3])){
544 digits->RemoveAt(i) ;
545 }
546 }
ba898748 547 }
454c62a4 548
ba898748 549 digits->Compress() ;
550 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
551 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
552 digit->SetIndexInList(i) ;
553 }
ba898748 554}
9a2cdbdf 555
d15a28e7 556//____________________________________________________________________________
9a1398dd 557Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
d15a28e7 558{
b2a60966 559 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
560
9f616d61 561 Bool_t rv = kFALSE ;
d15a28e7 562
9a2cdbdf 563 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
d15a28e7 564
2b629790 565 if(digit->GetId() <= nEMC ) rv = kTRUE;
ed4205d8 566
567 return rv ;
568}
569
ed4205d8 570//____________________________________________________________________________
9a1398dd 571Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
ed4205d8 572{
fad3e5b9 573 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
ed4205d8 574
575 Bool_t rv = kFALSE ;
88cb7938 576
9a2cdbdf 577 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
ed4205d8 578
2b629790 579 if(digit->GetId() > nEMC ) rv = kTRUE;
d15a28e7 580
581 return rv ;
582}
9a1398dd 583
88cb7938 584//____________________________________________________________________________
585void AliPHOSClusterizerv1::WriteRecPoints()
a4e98857 586{
7b7c1533 587
588 // Creates new branches with given title
a4e98857 589 // fills and writes into TreeR.
88cb7938 590
9a1398dd 591 Int_t index ;
ba898748 592 //Evaluate position, dispersion and other RecPoint properties..
9a2cdbdf 593 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
ba898748 594 for(index = 0; index < nEmc; index++){
cd79ec76 595 AliPHOSEmcRecPoint * rp =
9a2cdbdf 596 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
ba898748 597 rp->Purify(fEmcMinE) ;
a978bcea 598 if(rp->GetMultiplicity()==0){
9a2cdbdf 599 fEMCRecPoints->RemoveAt(index) ;
ba898748 600 delete rp ;
cd79ec76 601 continue;
ba898748 602 }
a978bcea 603
cd79ec76 604 // No vertex is available now, calculate corrections in PID
9a2cdbdf 605 rp->EvalAll(fW0,fDigitsArr) ;
cd79ec76 606 TVector3 fakeVtx(0.,0.,0.) ;
9a2cdbdf 607 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
cd79ec76 608 rp->EvalLocal2TrackingCSTransform();
ba898748 609 }
9a2cdbdf 610 fEMCRecPoints->Compress() ;
6d0a5c2e 611 fEMCRecPoints->Sort() ;
9a2cdbdf 612 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
613 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
614 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
ba898748 615 }
fbf811ec 616
af0b0570 617 //For each rec.point set the distance to the nearest bad crystal (BVP)
618 SetDistancesToBadChannels();
619
9a1398dd 620 //Now the same for CPV
9a2cdbdf 621 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
622 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
623 rp->EvalAll(fW0CPV,fDigitsArr) ;
76e42795 624 rp->EvalLocal2TrackingCSTransform();
ba898748 625 }
6d0a5c2e 626 fCPVRecPoints->Sort() ;
fbf811ec 627
9a2cdbdf 628 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
629 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
fbf811ec 630
9a2cdbdf 631 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
9a1398dd 632
ba898748 633 if(fWrite){ //We write TreeR
9a2cdbdf 634 fTreeR->Fill();
ba898748 635 }
9a1398dd 636}
637
638//____________________________________________________________________________
639void AliPHOSClusterizerv1::MakeClusters()
640{
641 // Steering method to construct the clusters stored in a list of Reconstructed Points
642 // A cluster is defined as a list of neighbour digits
d15a28e7 643
ba898748 644 //Remove digits below threshold
9a2cdbdf 645 CleanDigits(fDigitsArr) ;
ba898748 646
9a2cdbdf 647 TClonesArray * digitsC = static_cast<TClonesArray*>( fDigitsArr->Clone() ) ;
7b7c1533 648
d15a28e7 649 // Clusterization starts
9a1398dd 650
7b7c1533 651 TIter nextdigit(digitsC) ;
d15a28e7 652 AliPHOSDigit * digit ;
92862013 653 Bool_t notremoved = kTRUE ;
6ad0bfa0 654
88cb7938 655 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
656
657
9a1398dd 658 AliPHOSRecPoint * clu = 0 ;
c161df70 659
d1de15f5 660 TArrayI clusterdigitslist(1500) ;
d15a28e7 661 Int_t index ;
f2bc1b87 662
e95226ae 663 if (( IsInEmc (digit) &&
664 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
665 ( IsInCpv (digit) &&
666 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
d15a28e7 667 Int_t iDigitInCluster = 0 ;
7b7c1533 668
6ad0bfa0 669 if ( IsInEmc(digit) ) {
88cb7938 670 // start a new EMC RecPoint
9a2cdbdf 671 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
672 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
88cb7938 673
9a2cdbdf 674 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
675 clu = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
e95226ae 676 fNumberOfEmcClusters++ ;
677 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
678 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
679 iDigitInCluster++ ;
88cb7938 680 digitsC->Remove(digit) ;
ca77b032 681
d72dfbc3 682 } else {
88cb7938 683
684 // start a new CPV cluster
9a2cdbdf 685 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
686 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
88cb7938 687
9a2cdbdf 688 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
88cb7938 689
9a2cdbdf 690 clu = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
88cb7938 691 fNumberOfCpvClusters++ ;
e95226ae 692 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
88cb7938 693 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
694 iDigitInCluster++ ;
695 digitsC->Remove(digit) ;
d15a28e7 696 nextdigit.Reset() ;
88cb7938 697
698 // Here we remove remaining EMC digits, which cannot make a cluster
699
92862013 700 if( notremoved ) {
88cb7938 701 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
9f616d61 702 if( IsInEmc(digit) )
88cb7938 703 digitsC->Remove(digit) ;
d15a28e7 704 else
88cb7938 705 break ;
706 }
707 notremoved = kFALSE ;
708 }
709
d15a28e7 710 } // else
711
712 nextdigit.Reset() ;
fad3e5b9 713
d15a28e7 714 AliPHOSDigit * digitN ;
715 index = 0 ;
9f616d61 716 while (index < iDigitInCluster){ // scan over digits already in cluster
9a2cdbdf 717 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
88cb7938 718 index++ ;
719 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
720 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
d15a28e7 721 switch (ineb ) {
9f616d61 722 case 0 : // not a neighbour
88cb7938 723 break ;
724 case 1 : // are neighbours
e95226ae 725 if (IsInEmc (digitN))
726 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
727 else
728 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
729
88cb7938 730 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
731 iDigitInCluster++ ;
732 digitsC->Remove(digitN) ;
733 break ;
9f616d61 734 case 2 : // too far from each other
88cb7938 735 goto endofloop;
736 } // switch
737
738 } // while digitN
739
d15a28e7 740 endofloop: ;
88cb7938 741 nextdigit.Reset() ;
742
d15a28e7 743 } // loop over cluster
ad8cfaf4 744
ad8cfaf4 745 } // energy theshold
ad8cfaf4 746
31aa6d6c 747
d15a28e7 748 } // while digit
749
7b7c1533 750 delete digitsC ;
9688c1dd 751
9a1398dd 752}
753
754//____________________________________________________________________________
a4e98857 755void AliPHOSClusterizerv1::MakeUnfolding()
756{
757 // Unfolds clusters using the shape of an ElectroMagnetic shower
9688c1dd 758 // Performs unfolding of all EMC/CPV clusters
9a1398dd 759
a4e98857 760 // Unfold first EMC clusters
9a1398dd 761 if(fNumberOfEmcClusters > 0){
762
9a2cdbdf 763 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 764
765 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
766 Int_t index ;
767 for(index = 0 ; index < numberofNotUnfolded ; index++){
768
9a2cdbdf 769 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
9a1398dd 770 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 771 break ;
9a1398dd 772
773 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 774 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 775 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 776 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
9a1398dd 777
778 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 779 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
9a2cdbdf 780 fEMCRecPoints->Remove(emcRecPoint);
781 fEMCRecPoints->Compress() ;
88cb7938 782 index-- ;
783 fNumberOfEmcClusters -- ;
784 numberofNotUnfolded-- ;
9a1398dd 785 }
c6bf27f2 786 else{
787 emcRecPoint->SetNExMax(1) ; //Only one local maximum
788 }
9a1398dd 789
790 delete[] maxAt ;
791 delete[] maxAtEnergy ;
792 }
793 }
a4e98857 794 // Unfolding of EMC clusters finished
9a1398dd 795
796
a4e98857 797 // Unfold now CPV clusters
9a1398dd 798 if(fNumberOfCpvClusters > 0){
799
9a2cdbdf 800 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 801
802 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
803 Int_t index ;
804 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
805
9a2cdbdf 806 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
9a1398dd 807
808 if(recPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 809 break ;
9a1398dd 810
88cb7938 811 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
9a1398dd 812
813 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 814 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 815 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 816 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
9a1398dd 817
818 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 819 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
9a2cdbdf 820 fCPVRecPoints->Remove(emcRecPoint);
821 fCPVRecPoints->Compress() ;
88cb7938 822 index-- ;
823 numberofCpvNotUnfolded-- ;
824 fNumberOfCpvClusters-- ;
9a1398dd 825 }
826
827 delete[] maxAt ;
828 delete[] maxAtEnergy ;
829 }
830 }
831 //Unfolding of Cpv clusters finished
832
833}
834
9a1398dd 835//____________________________________________________________________________
a978bcea 836Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
9a1398dd 837{
838 // Shape of the shower (see PHOS TDR)
a4e98857 839 // If you change this function, change also the gradient evaluation in ChiSquare()
9a1398dd 840
a978bcea 841 //for the moment we neglect dependence on the incident angle.
842
843 Double_t r2 = x*x + z*z ;
844 Double_t r4 = r2*r2 ;
845 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
9a1398dd 846 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
847 return shape ;
848}
849
850//____________________________________________________________________________
851void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
88cb7938 852 Int_t nMax,
853 AliPHOSDigit ** maxAt,
854 Float_t * maxAtEnergy)
9a1398dd 855{
856 // Performs the unfolding of a cluster with nMax overlapping showers
857
858 Int_t nPar = 3 * nMax ;
859 Float_t * fitparameters = new Float_t[nPar] ;
860
861 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
862 if( !rv ) {
88cb7938 863 // Fit failed, return and remove cluster
c6bf27f2 864 iniEmc->SetNExMax(-1) ;
9a1398dd 865 delete[] fitparameters ;
866 return ;
867 }
868
869 // create ufolded rec points and fill them with new energy lists
870 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
871 // and later correct this number in acordance with actual energy deposition
872
873 Int_t nDigits = iniEmc->GetMultiplicity() ;
874 Float_t * efit = new Float_t[nDigits] ;
a978bcea 875 Float_t xDigit=0.,zDigit=0. ;
7b7c1533 876 Float_t xpar=0.,zpar=0.,epar=0. ;
9a1398dd 877 Int_t relid[4] ;
7b7c1533 878 AliPHOSDigit * digit = 0 ;
9a1398dd 879 Int_t * emcDigits = iniEmc->GetDigitsList() ;
880
a978bcea 881 TVector3 vIncid ;
882
9a1398dd 883 Int_t iparam ;
884 Int_t iDigit ;
885 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 886 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
887 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
888 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
9a1398dd 889 efit[iDigit] = 0;
890
891 iparam = 0 ;
892 while(iparam < nPar ){
893 xpar = fitparameters[iparam] ;
894 zpar = fitparameters[iparam+1] ;
895 epar = fitparameters[iparam+2] ;
896 iparam += 3 ;
9a2cdbdf 897// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
a978bcea 898// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
899 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 900 }
901 }
902
903
904 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
905 // so that energy deposited in each cell is distributed betwin new clusters proportionally
906 // to its contribution to efit
907
908 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
909 Float_t ratio ;
910
911 iparam = 0 ;
912 while(iparam < nPar ){
913 xpar = fitparameters[iparam] ;
914 zpar = fitparameters[iparam+1] ;
915 epar = fitparameters[iparam+2] ;
916 iparam += 3 ;
9a2cdbdf 917// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
9a1398dd 918
7b7c1533 919 AliPHOSEmcRecPoint * emcRP = 0 ;
9a1398dd 920
9a2cdbdf 921 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
9a1398dd 922
9a2cdbdf 923 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
924 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
9a1398dd 925
9a2cdbdf 926 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
927 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
9a1398dd 928 fNumberOfEmcClusters++ ;
c6bf27f2 929 emcRP->SetNExMax((Int_t)nPar/3) ;
9a1398dd 930 }
9a2cdbdf 931 else{//create new entries in fCPVRecPoints
932 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
933 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
9a1398dd 934
9a2cdbdf 935 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
936 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
9a1398dd 937 fNumberOfCpvClusters++ ;
938 }
939
940 Float_t eDigit ;
941 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 942 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
943 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
944 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
a978bcea 945// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
946 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
9a1398dd 947 eDigit = emcEnergies[iDigit] * ratio ;
948 emcRP->AddDigit( *digit, eDigit ) ;
88cb7938 949 }
9a1398dd 950 }
951
952 delete[] fitparameters ;
953 delete[] efit ;
954
955}
956
957//_____________________________________________________________________________
958void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
959{
a4e98857 960 // Calculates the Chi square for the cluster unfolding minimization
9a1398dd 961 // Number of parameters, Gradient, Chi squared, parameters, what to do
962
88cb7938 963 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
9a1398dd 964
88cb7938 965 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
966 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
9a2cdbdf 967 // A bit buggy way to get an access to the geometry
968 // To be revised!
969 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
970
971// TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
9a1398dd 972
88cb7938 973 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
9a1398dd 974
975 Int_t * emcDigits = emcRP->GetDigitsList() ;
976
7b7c1533 977 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
9a1398dd 978
979 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
88cb7938 980
a978bcea 981// TVector3 vInc ;
9a1398dd 982 fret = 0. ;
983 Int_t iparam ;
984
985 if(iflag == 2)
986 for(iparam = 0 ; iparam < nPar ; iparam++)
987 Grad[iparam] = 0 ; // Will evaluate gradient
988
989 Double_t efit ;
990
991 AliPHOSDigit * digit ;
992 Int_t iDigit ;
993
7b7c1533 994 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
9a1398dd 995
88cb7938 996 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
9a1398dd 997
998 Int_t relid[4] ;
999 Float_t xDigit ;
1000 Float_t zDigit ;
1001
1002 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1003
1004 geom->RelPosInModule(relid, xDigit, zDigit) ;
1005
1006 if(iflag == 2){ // calculate gradient
1007 Int_t iParam = 0 ;
1008 efit = 0 ;
1009 while(iParam < nPar ){
a978bcea 1010 Double_t dx = (xDigit - x[iParam]) ;
88cb7938 1011 iParam++ ;
a978bcea 1012 Double_t dz = (zDigit - x[iParam]) ;
88cb7938 1013 iParam++ ;
9a2cdbdf 1014// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
a978bcea 1015// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
1016 efit += x[iParam] * ShowerShape(dx,dz) ;
88cb7938 1017 iParam++ ;
9a1398dd 1018 }
88cb7938 1019 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
9a1398dd 1020 iParam = 0 ;
1021 while(iParam < nPar ){
88cb7938 1022 Double_t xpar = x[iParam] ;
1023 Double_t zpar = x[iParam+1] ;
1024 Double_t epar = x[iParam+2] ;
9a2cdbdf 1025// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
88cb7938 1026 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
a978bcea 1027// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1028 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1029//DP: No incident angle dependence in gradient yet!!!!!!
88cb7938 1030 Double_t r4 = dr*dr*dr*dr ;
1031 Double_t r295 = TMath::Power(dr,2.95) ;
1032 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1033 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1034
1035 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1036 iParam++ ;
1037 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1038 iParam++ ;
1039 Grad[iParam] += shape ; // Derivative over energy
1040 iParam++ ;
9a1398dd 1041 }
1042 }
1043 efit = 0;
1044 iparam = 0 ;
1045
1046 while(iparam < nPar ){
1047 Double_t xpar = x[iparam] ;
1048 Double_t zpar = x[iparam+1] ;
1049 Double_t epar = x[iparam+2] ;
1050 iparam += 3 ;
9a2cdbdf 1051// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
a978bcea 1052// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1053 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 1054 }
1055
88cb7938 1056 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
9a1398dd 1057 // Here we assume, that sigma = sqrt(E)
1058 }
83974468 1059
d15a28e7 1060}
1061
1062//____________________________________________________________________________
702ab87e 1063void AliPHOSClusterizerv1::Print(const Option_t *)const
d15a28e7 1064{
d1de15f5 1065 // Print clusterizer parameters
1066
21cd0c07 1067 TString message ;
1068 TString taskName(GetName()) ;
1069 taskName.ReplaceAll(Version(), "") ;
1070
1071 if( strcmp(GetName(), "") !=0 ) {
9a1398dd 1072 // Print parameters
21cd0c07 1073 message = "\n--------------- %s %s -----------\n" ;
1074 message += "Clusterizing digits from the file: %s\n" ;
1075 message += " Branch: %s\n" ;
1076 message += " EMC Clustering threshold = %f\n" ;
1077 message += " EMC Local Maximum cut = %f\n" ;
1078 message += " EMC Logarothmic weight = %f\n" ;
1079 message += " CPV Clustering threshold = %f\n" ;
1080 message += " CPV Local Maximum cut = %f\n" ;
1081 message += " CPV Logarothmic weight = %f\n" ;
9a1398dd 1082 if(fToUnfold)
21cd0c07 1083 message += " Unfolding on\n" ;
9a1398dd 1084 else
21cd0c07 1085 message += " Unfolding off\n" ;
9a1398dd 1086
21cd0c07 1087 message += "------------------------------------------------------------------" ;
9a1398dd 1088 }
21cd0c07 1089 else
1090 message = " AliPHOSClusterizerv1 not initialized " ;
1091
351dd634 1092 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
21cd0c07 1093 taskName.Data(),
1094 GetTitle(),
1095 taskName.Data(),
1096 GetName(),
1097 fEmcClusteringThreshold,
1098 fEmcLocMaxCut,
1099 fW0,
1100 fCpvClusteringThreshold,
1101 fCpvLocMaxCut,
351dd634 1102 fW0CPV )) ;
9a1398dd 1103}
a978bcea 1104//____________________________________________________________________________
a4e98857 1105void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1106{
1107 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
9a1398dd 1108
9a2cdbdf 1109 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1110 fEMCRecPoints->GetEntriesFast(),
1111 fCPVRecPoints->GetEntriesFast() )) ;
88cb7938 1112
9a1398dd 1113 if(strstr(option,"all")) {
709e117a 1114 printf("\n EMC clusters \n") ;
1115 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
9a1398dd 1116 Int_t index ;
9a2cdbdf 1117 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1118 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
9a1398dd 1119 TVector3 locpos;
1120 rp->GetLocalPosition(locpos);
9a1398dd 1121 Float_t lambda[2];
1122 rp->GetElipsAxis(lambda);
9a1398dd 1123 Int_t * primaries;
1124 Int_t nprimaries;
1125 primaries = rp->GetPrimaries(nprimaries);
709e117a 1126 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
11f9c5ff 1127 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1128 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
3bf72d32 1129
21cd0c07 1130 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
709e117a 1131 printf("%d ", primaries[iprimary] ) ;
21cd0c07 1132 }
709e117a 1133 printf("\n") ;
1134 }
1135
1136 //Now plot CPV recPoints
1137 printf("\n CPV clusters \n") ;
1138 printf("Index Ene(MeV) Module X Y Z \n") ;
9a2cdbdf 1139 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1140 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
2bb500e5 1141
1142 TVector3 locpos;
1143 rp->GetLocalPosition(locpos);
1144
1145 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1146 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1147 locpos.X(), locpos.Y(), locpos.Z()) ;
88cb7938 1148 }
9a1398dd 1149 }
d15a28e7 1150}
9a1398dd 1151
af0b0570 1152
1153//____________________________________________________________________________
1154void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1155{
1156 //For each EMC rec. point set the distance to the nearest bad crystal.
1157 //Author: Boris Polichtchouk
1158
e68222ce 1159 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1160 AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
af0b0570 1161
af0b0570 1162 Int_t badIds[8000];
e68222ce 1163 fgCalibData->EmcBadChannelIds(badIds);
af0b0570 1164
1165 AliPHOSEmcRecPoint* rp;
1166
1167 TMatrixF gmat;
1168 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1169 TVector3 gposBadChannel; // global position of bad crystal
1170 TVector3 dR;
1171
1172 Float_t dist,minDist;
1173
9a2cdbdf 1174 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1175 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
af0b0570 1176 minDist = 1.e+07;
1177
e68222ce 1178 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
af0b0570 1179 rp->GetGlobalPosition(gposRecPoint,gmat);
9a2cdbdf 1180 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
af0b0570 1181 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1182 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1183 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1184 dR = gposBadChannel-gposRecPoint;
1185 dist = dR.Mag();
1186 if(dist<minDist) minDist = dist;
1187 }
1188
1189 rp->SetDistanceToBadCrystal(minDist);
1190 }
1191
1192}