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