AliEveEventManager::AssertGeometry() - consolidate with changes in
[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),
d07a44a8 197 fWrite(0),
198 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
a28333c5 199 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
0378398c 200 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
a28333c5 201 fW0CPV(0), fEmcTimeGate(0)
d15a28e7 202{
f035f6ce 203 // default ctor (to be used mainly by Streamer)
f035f6ce 204
8d0f3f77 205 InitParameters() ;
92f521a9 206 fDefaultInit = kTRUE ;
9a1398dd 207}
7b7c1533 208
9a1398dd 209//____________________________________________________________________________
9a2cdbdf 210AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
211 AliPHOSClusterizer(geom),
0378398c 212 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
d07a44a8 213 fWrite(0),
214 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
a28333c5 215 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
0378398c 216 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
a28333c5 217 fW0CPV(0), fEmcTimeGate(0)
9a1398dd 218{
a4e98857 219 // ctor with the indication of the file where header Tree and digits Tree are stored
9a1398dd 220
8d0f3f77 221 InitParameters() ;
2bd5457f 222 Init() ;
92f521a9 223 fDefaultInit = kFALSE ;
9a1398dd 224}
fc12304f 225
9688c1dd 226//____________________________________________________________________________
227 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
228{
0bc3b8ed 229 // dtor
230
9688c1dd 231}
3758d9fc 232//____________________________________________________________________________
9a2cdbdf 233void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
a4e98857 234{
9a2cdbdf 235 // Steering method to perform clusterization for one event
236 // The input is the tree with digits.
237 // The output is the tree with clusters.
9a1398dd 238
239 if(strstr(option,"tim"))
7d493c2c 240 gBenchmark->Start("PHOSClusterizer");
9a1398dd 241
eabde521 242 if(strstr(option,"print")) {
88cb7938 243 Print() ;
eabde521 244 return ;
245 }
bed9e3fb 246
9a2cdbdf 247 MakeClusters() ;
9a1398dd 248
9a2cdbdf 249 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
250 fEMCRecPoints->GetEntries()));
251 if(AliLog::GetGlobalDebugLevel()>1)
252 fEMCRecPoints->Print();
88cb7938 253
9a2cdbdf 254 if(fToUnfold)
255 MakeUnfolding();
ddd1a39c 256
9a2cdbdf 257 WriteRecPoints();
7b7c1533 258
9a2cdbdf 259 if(strstr(option,"deb"))
260 PrintRecPoints(option) ;
94de8339 261
9a1398dd 262 if(strstr(option,"tim")){
263 gBenchmark->Stop("PHOSClusterizer");
9a2cdbdf 264 AliInfo(Form("took %f seconds for Clusterizing\n",
265 gBenchmark->GetCpuTime("PHOSClusterizer")));
8d8258f6 266 }
9a1398dd 267}
268
269//____________________________________________________________________________
a0636361 270Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
d1de15f5 271 Int_t nPar, Float_t * fitparameters) const
9a1398dd 272{
273 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
a4e98857 274 // The initial values for fitting procedure are set equal to the positions of local maxima.
f035f6ce 275 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
9a1398dd 276
88cb7938 277
9a1398dd 278 gMinuit->mncler(); // Reset Minuit's list of paramters
279 gMinuit->SetPrintLevel(-1) ; // No Printout
280 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
281 // To set the address of the minimization function
282
283 TList * toMinuit = new TList();
284 toMinuit->AddAt(emcRP,0) ;
9a2cdbdf 285 toMinuit->AddAt(fDigitsArr,1) ;
286 toMinuit->AddAt(fGeom,2) ;
9a1398dd 287
288 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
289
290 // filling initial values for fit parameters
291 AliPHOSDigit * digit ;
292
293 Int_t ierflg = 0;
294 Int_t index = 0 ;
295 Int_t nDigits = (Int_t) nPar / 3 ;
296
297 Int_t iDigit ;
298
9a1398dd 299 for(iDigit = 0; iDigit < nDigits; iDigit++){
a0636361 300 digit = maxAt[iDigit];
9a1398dd 301
302 Int_t relid[4] ;
7b7c1533 303 Float_t x = 0.;
304 Float_t z = 0.;
9a2cdbdf 305 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
306 fGeom->RelPosInModule(relid, x, z) ;
9a1398dd 307
308 Float_t energy = maxAtEnergy[iDigit] ;
309
310 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
311 index++ ;
312 if(ierflg != 0){
21cd0c07 313 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
9a1398dd 314 return kFALSE;
315 }
316 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
317 index++ ;
318 if(ierflg != 0){
21cd0c07 319 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
9a1398dd 320 return kFALSE;
321 }
322 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
323 index++ ;
324 if(ierflg != 0){
21cd0c07 325 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
9a1398dd 326 return kFALSE;
327 }
328 }
329
330 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
331 // depends on it.
332 Double_t p1 = 1.0 ;
333 Double_t p2 = 0.0 ;
334
335 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
336 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
337 gMinuit->SetMaxIterations(5);
338 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
339
340 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
341
342 if(ierflg == 4){ // Minimum not found
21cd0c07 343 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
9a1398dd 344 return kFALSE ;
345 }
346 for(index = 0; index < nPar; index++){
347 Double_t err ;
348 Double_t val ;
349 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
350 fitparameters[index] = val ;
351 }
352
353 delete toMinuit ;
354 return kTRUE;
c3c187e5 355
d15a28e7 356}
357
548f0134 358
d15a28e7 359//____________________________________________________________________________
a4e98857 360void AliPHOSClusterizerv1::Init()
361{
2bd5457f 362 // Make all memory allocations which can not be done in default constructor.
363 // Attach the Clusterizer task to the list of PHOS tasks
ba898748 364
9a2cdbdf 365 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
3758d9fc 366
7b7c1533 367 if(!gMinuit)
88cb7938 368 gMinuit = new TMinuit(100);
fbf811ec 369
a28333c5 370 if (!fgCalibData)
371 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
372 if (fgCalibData->GetCalibDataEmc() == 0)
373 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
e68222ce 374
9a1398dd 375}
7b7c1533 376
8d0f3f77 377//____________________________________________________________________________
378void AliPHOSClusterizerv1::InitParameters()
379{
380
381 fNumberOfCpvClusters = 0 ;
382 fNumberOfEmcClusters = 0 ;
3799bcb5 383
384 const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
385 if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
386
387 const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv();
388 if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
389
390 fCpvClusteringThreshold = parCpv->GetClusteringThreshold();
391 fEmcClusteringThreshold = parEmc->GetClusteringThreshold();
8d0f3f77 392
3799bcb5 393 fEmcLocMaxCut = parEmc->GetLocalMaxCut();
394 fCpvLocMaxCut = parCpv->GetLocalMaxCut();
ba898748 395
3799bcb5 396 fW0 = parEmc->GetLogWeight();
397 fW0CPV = parCpv->GetLogWeight();
8d0f3f77 398
454c62a4 399 fEmcTimeGate = 1.e-6 ;
8d0f3f77 400
a42611a9 401 fToUnfold = parEmc->ToUnfold() ;
88cb7938 402
ba898748 403 fWrite = kTRUE ;
8d0f3f77 404}
405
9a1398dd 406//____________________________________________________________________________
407Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
d15a28e7 408{
b2a60966 409 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
410 // = 1 are neighbour
411 // = 2 are not neighbour but do not continue searching
d07a44a8 412 // =-1 are not neighbour, continue searching, but do not look before d2 next time
a4e98857 413 // neighbours are defined as digits having at least a common vertex
b2a60966 414 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
415 // which is compared to a digit (d2) not yet in a cluster
416
d15a28e7 417 Int_t relid1[4] ;
9a2cdbdf 418 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
d15a28e7 419
420 Int_t relid2[4] ;
9a2cdbdf 421 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
d15a28e7 422
9688c1dd 423 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
92862013 424 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
425 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
d15a28e7 426
454c62a4 427 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
428 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
9688c1dd 429 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
d07a44a8 430 return 1 ;
d15a28e7 431 }
432 else {
433 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
d07a44a8 434 return 2; // Difference in row numbers is too large to look further
d15a28e7 435 }
d07a44a8 436 return 0 ;
d15a28e7 437
438 }
439 else {
d07a44a8 440 if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
441 return -1 ;
442 if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
443 return -1 ;
d15a28e7 444
d07a44a8 445 return 2 ;
d15a28e7 446
447 }
d72dfbc3 448
d07a44a8 449 return 0 ;
d15a28e7 450}
d15a28e7 451//____________________________________________________________________________
9a1398dd 452Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
d15a28e7 453{
b2a60966 454 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
455
9f616d61 456 Bool_t rv = kFALSE ;
d15a28e7 457
9a2cdbdf 458 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
d15a28e7 459
2b629790 460 if(digit->GetId() <= nEMC ) rv = kTRUE;
ed4205d8 461
462 return rv ;
463}
464
ed4205d8 465//____________________________________________________________________________
9a1398dd 466Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
ed4205d8 467{
fad3e5b9 468 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
ed4205d8 469
470 Bool_t rv = kFALSE ;
88cb7938 471
9a2cdbdf 472 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
ed4205d8 473
2b629790 474 if(digit->GetId() > nEMC ) rv = kTRUE;
d15a28e7 475
476 return rv ;
477}
9a1398dd 478
88cb7938 479//____________________________________________________________________________
480void AliPHOSClusterizerv1::WriteRecPoints()
a4e98857 481{
7b7c1533 482
483 // Creates new branches with given title
a4e98857 484 // fills and writes into TreeR.
88cb7938 485
9a1398dd 486 Int_t index ;
ba898748 487 //Evaluate position, dispersion and other RecPoint properties..
9a2cdbdf 488 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
a28333c5 489 Float_t emcMinE= AliPHOSReconstructor::GetRecoParamEmc()->GetMinE(); //Minimal digit energy
ba898748 490 for(index = 0; index < nEmc; index++){
cd79ec76 491 AliPHOSEmcRecPoint * rp =
9a2cdbdf 492 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
a28333c5 493 rp->Purify(emcMinE) ;
a978bcea 494 if(rp->GetMultiplicity()==0){
9a2cdbdf 495 fEMCRecPoints->RemoveAt(index) ;
ba898748 496 delete rp ;
cd79ec76 497 continue;
ba898748 498 }
a978bcea 499
cd79ec76 500 // No vertex is available now, calculate corrections in PID
9a2cdbdf 501 rp->EvalAll(fW0,fDigitsArr) ;
cd79ec76 502 TVector3 fakeVtx(0.,0.,0.) ;
9a2cdbdf 503 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
cd79ec76 504 rp->EvalLocal2TrackingCSTransform();
ba898748 505 }
9a2cdbdf 506 fEMCRecPoints->Compress() ;
6d0a5c2e 507 fEMCRecPoints->Sort() ;
9a2cdbdf 508 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
509 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
510 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
ba898748 511 }
fbf811ec 512
af0b0570 513 //For each rec.point set the distance to the nearest bad crystal (BVP)
514 SetDistancesToBadChannels();
515
9a1398dd 516 //Now the same for CPV
9a2cdbdf 517 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
518 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
519 rp->EvalAll(fW0CPV,fDigitsArr) ;
76e42795 520 rp->EvalLocal2TrackingCSTransform();
ba898748 521 }
6d0a5c2e 522 fCPVRecPoints->Sort() ;
fbf811ec 523
9a2cdbdf 524 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
525 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
fbf811ec 526
9a2cdbdf 527 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
9a1398dd 528
ba898748 529 if(fWrite){ //We write TreeR
9a2cdbdf 530 fTreeR->Fill();
ba898748 531 }
9a1398dd 532}
533
534//____________________________________________________________________________
535void AliPHOSClusterizerv1::MakeClusters()
536{
537 // Steering method to construct the clusters stored in a list of Reconstructed Points
538 // A cluster is defined as a list of neighbour digits
d15a28e7 539
8d8258f6 540 fNumberOfCpvClusters = 0 ;
541 fNumberOfEmcClusters = 0 ;
542
d07a44a8 543 //Mark all digits as unused yet
544 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
545 for(Int_t i=0; i<nDigits; i++){
546 fDigitsUsed[i]=0 ;
547 }
548 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
549 //e.g. first digit in this module, first CPV digit etc.
d15a28e7 550 AliPHOSDigit * digit ;
d07a44a8 551 TArrayI clusterdigitslist(1500) ;
552 AliPHOSRecPoint * clu = 0 ;
553 for(Int_t i=0; i<nDigits; i++){
554 if(fDigitsUsed[i])
555 continue ;
6ad0bfa0 556
d07a44a8 557 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
88cb7938 558
d07a44a8 559 clu=0 ;
c161df70 560
d15a28e7 561 Int_t index ;
f2bc1b87 562
d07a44a8 563 //is this digit so energetic that start cluster?
a28333c5 564 if (( IsInEmc(digit) && digit->GetEnergy() > fEmcClusteringThreshold ) ||
565 ( IsInCpv(digit) && digit->GetEnergy() > fCpvClusteringThreshold ) ) {
d15a28e7 566 Int_t iDigitInCluster = 0 ;
7b7c1533 567
6ad0bfa0 568 if ( IsInEmc(digit) ) {
88cb7938 569 // start a new EMC RecPoint
9a2cdbdf 570 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
571 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
88cb7938 572
9a2cdbdf 573 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
d07a44a8 574 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
e95226ae 575 fNumberOfEmcClusters++ ;
a28333c5 576 clu->AddDigit(*digit, digit->GetEnergy()) ;
e95226ae 577 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
578 iDigitInCluster++ ;
d07a44a8 579 fDigitsUsed[i]=kTRUE ;
d72dfbc3 580 } else {
88cb7938 581 // start a new CPV cluster
9a2cdbdf 582 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
583 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
88cb7938 584
9a2cdbdf 585 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
d07a44a8 586 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
88cb7938 587 fNumberOfCpvClusters++ ;
a28333c5 588 clu->AddDigit(*digit, digit->GetEnergy()) ;
88cb7938 589 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
590 iDigitInCluster++ ;
d07a44a8 591 fDigitsUsed[i]=kTRUE ;
d15a28e7 592 } // else
593
d07a44a8 594 //Now scan remaining digits in list to find neigbours of our seed
595
d15a28e7 596 AliPHOSDigit * digitN ;
597 index = 0 ;
9f616d61 598 while (index < iDigitInCluster){ // scan over digits already in cluster
d07a44a8 599 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
88cb7938 600 index++ ;
d07a44a8 601 for(Int_t j=iFirst; j<nDigits; j++){
602 if(fDigitsUsed[j])
603 continue ; //look through remaining digits
604 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
88cb7938 605 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
d15a28e7 606 switch (ineb ) {
d07a44a8 607 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
608 iFirst=j ;
609 break ;
9f616d61 610 case 0 : // not a neighbour
88cb7938 611 break ;
612 case 1 : // are neighbours
a28333c5 613 clu->AddDigit(*digitN, digitN->GetEnergy());
d07a44a8 614 clusterdigitslist[iDigitInCluster] = j ;
88cb7938 615 iDigitInCluster++ ;
d07a44a8 616 fDigitsUsed[j]=kTRUE ;
88cb7938 617 break ;
9f616d61 618 case 2 : // too far from each other
88cb7938 619 goto endofloop;
620 } // switch
621
d07a44a8 622 }
88cb7938 623
d07a44a8 624 endofloop: ; //scanned all possible neighbours for this digit
88cb7938 625
d15a28e7 626 } // loop over cluster
ad8cfaf4 627 } // energy theshold
d07a44a8 628 }
9688c1dd 629
9a1398dd 630}
631
632//____________________________________________________________________________
a4e98857 633void AliPHOSClusterizerv1::MakeUnfolding()
634{
635 // Unfolds clusters using the shape of an ElectroMagnetic shower
9688c1dd 636 // Performs unfolding of all EMC/CPV clusters
9a1398dd 637
a4e98857 638 // Unfold first EMC clusters
9a1398dd 639 if(fNumberOfEmcClusters > 0){
640
9a2cdbdf 641 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 642
643 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
644 Int_t index ;
645 for(index = 0 ; index < numberofNotUnfolded ; index++){
646
9a2cdbdf 647 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
9a1398dd 648 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 649 break ;
9a1398dd 650
651 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 652 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 653 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 654 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
9a1398dd 655
656 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 657 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
a28333c5 658
9a2cdbdf 659 fEMCRecPoints->Remove(emcRecPoint);
660 fEMCRecPoints->Compress() ;
88cb7938 661 index-- ;
662 fNumberOfEmcClusters -- ;
663 numberofNotUnfolded-- ;
9a1398dd 664 }
c6bf27f2 665 else{
666 emcRecPoint->SetNExMax(1) ; //Only one local maximum
667 }
9a1398dd 668
669 delete[] maxAt ;
670 delete[] maxAtEnergy ;
671 }
672 }
a4e98857 673 // Unfolding of EMC clusters finished
9a1398dd 674
675
a4e98857 676 // Unfold now CPV clusters
9a1398dd 677 if(fNumberOfCpvClusters > 0){
678
9a2cdbdf 679 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 680
681 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
682 Int_t index ;
683 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
684
9a2cdbdf 685 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
9a1398dd 686
687 if(recPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 688 break ;
9a1398dd 689
88cb7938 690 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
9a1398dd 691
692 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 693 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 694 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 695 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
9a1398dd 696
697 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 698 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
9a2cdbdf 699 fCPVRecPoints->Remove(emcRecPoint);
700 fCPVRecPoints->Compress() ;
88cb7938 701 index-- ;
702 numberofCpvNotUnfolded-- ;
703 fNumberOfCpvClusters-- ;
9a1398dd 704 }
705
706 delete[] maxAt ;
707 delete[] maxAtEnergy ;
708 }
709 }
710 //Unfolding of Cpv clusters finished
711
712}
713
9a1398dd 714//____________________________________________________________________________
a978bcea 715Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
9a1398dd 716{
717 // Shape of the shower (see PHOS TDR)
a4e98857 718 // If you change this function, change also the gradient evaluation in ChiSquare()
9a1398dd 719
a978bcea 720 //for the moment we neglect dependence on the incident angle.
721
722 Double_t r2 = x*x + z*z ;
723 Double_t r4 = r2*r2 ;
724 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
9a1398dd 725 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
726 return shape ;
727}
728
729//____________________________________________________________________________
730void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
88cb7938 731 Int_t nMax,
732 AliPHOSDigit ** maxAt,
733 Float_t * maxAtEnergy)
9a1398dd 734{
735 // Performs the unfolding of a cluster with nMax overlapping showers
736
737 Int_t nPar = 3 * nMax ;
738 Float_t * fitparameters = new Float_t[nPar] ;
739
740 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
a28333c5 741
9a1398dd 742 if( !rv ) {
88cb7938 743 // Fit failed, return and remove cluster
c6bf27f2 744 iniEmc->SetNExMax(-1) ;
9a1398dd 745 delete[] fitparameters ;
746 return ;
747 }
748
749 // create ufolded rec points and fill them with new energy lists
750 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
751 // and later correct this number in acordance with actual energy deposition
752
753 Int_t nDigits = iniEmc->GetMultiplicity() ;
754 Float_t * efit = new Float_t[nDigits] ;
a978bcea 755 Float_t xDigit=0.,zDigit=0. ;
7b7c1533 756 Float_t xpar=0.,zpar=0.,epar=0. ;
9a1398dd 757 Int_t relid[4] ;
7b7c1533 758 AliPHOSDigit * digit = 0 ;
9a1398dd 759 Int_t * emcDigits = iniEmc->GetDigitsList() ;
760
a978bcea 761 TVector3 vIncid ;
762
9a1398dd 763 Int_t iparam ;
764 Int_t iDigit ;
765 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 766 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
767 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
768 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
9a1398dd 769 efit[iDigit] = 0;
770
771 iparam = 0 ;
772 while(iparam < nPar ){
773 xpar = fitparameters[iparam] ;
774 zpar = fitparameters[iparam+1] ;
775 epar = fitparameters[iparam+2] ;
776 iparam += 3 ;
9a2cdbdf 777// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
a978bcea 778// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
779 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 780 }
781 }
782
9a1398dd 783 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
784 // so that energy deposited in each cell is distributed betwin new clusters proportionally
785 // to its contribution to efit
786
787 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
788 Float_t ratio ;
789
790 iparam = 0 ;
791 while(iparam < nPar ){
792 xpar = fitparameters[iparam] ;
793 zpar = fitparameters[iparam+1] ;
794 epar = fitparameters[iparam+2] ;
795 iparam += 3 ;
9a2cdbdf 796// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
a28333c5 797
7b7c1533 798 AliPHOSEmcRecPoint * emcRP = 0 ;
9a1398dd 799
9a2cdbdf 800 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
9a1398dd 801
9a2cdbdf 802 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
803 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
9a1398dd 804
9a2cdbdf 805 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
806 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
9a1398dd 807 fNumberOfEmcClusters++ ;
c6bf27f2 808 emcRP->SetNExMax((Int_t)nPar/3) ;
9a1398dd 809 }
9a2cdbdf 810 else{//create new entries in fCPVRecPoints
811 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
812 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
9a1398dd 813
9a2cdbdf 814 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
815 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
9a1398dd 816 fNumberOfCpvClusters++ ;
817 }
818
819 Float_t eDigit ;
820 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 821 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
822 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
823 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
a978bcea 824// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
825 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
9a1398dd 826 eDigit = emcEnergies[iDigit] * ratio ;
827 emcRP->AddDigit( *digit, eDigit ) ;
88cb7938 828 }
9a1398dd 829 }
830
831 delete[] fitparameters ;
832 delete[] efit ;
833
834}
835
836//_____________________________________________________________________________
837void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
838{
a4e98857 839 // Calculates the Chi square for the cluster unfolding minimization
9a1398dd 840 // Number of parameters, Gradient, Chi squared, parameters, what to do
841
88cb7938 842 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
9a1398dd 843
88cb7938 844 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
845 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
9a2cdbdf 846 // A bit buggy way to get an access to the geometry
847 // To be revised!
848 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
849
850// TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
9a1398dd 851
88cb7938 852 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
9a1398dd 853
854 Int_t * emcDigits = emcRP->GetDigitsList() ;
855
7b7c1533 856 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
9a1398dd 857
858 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
88cb7938 859
a978bcea 860// TVector3 vInc ;
9a1398dd 861 fret = 0. ;
862 Int_t iparam ;
863
864 if(iflag == 2)
865 for(iparam = 0 ; iparam < nPar ; iparam++)
866 Grad[iparam] = 0 ; // Will evaluate gradient
867
868 Double_t efit ;
869
870 AliPHOSDigit * digit ;
871 Int_t iDigit ;
872
7b7c1533 873 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
9a1398dd 874
88cb7938 875 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
9a1398dd 876
877 Int_t relid[4] ;
878 Float_t xDigit ;
879 Float_t zDigit ;
880
881 geom->AbsToRelNumbering(digit->GetId(), relid) ;
882
883 geom->RelPosInModule(relid, xDigit, zDigit) ;
884
885 if(iflag == 2){ // calculate gradient
886 Int_t iParam = 0 ;
887 efit = 0 ;
888 while(iParam < nPar ){
a978bcea 889 Double_t dx = (xDigit - x[iParam]) ;
88cb7938 890 iParam++ ;
a978bcea 891 Double_t dz = (zDigit - x[iParam]) ;
88cb7938 892 iParam++ ;
9a2cdbdf 893// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
a978bcea 894// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
895 efit += x[iParam] * ShowerShape(dx,dz) ;
88cb7938 896 iParam++ ;
9a1398dd 897 }
88cb7938 898 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
9a1398dd 899 iParam = 0 ;
900 while(iParam < nPar ){
88cb7938 901 Double_t xpar = x[iParam] ;
902 Double_t zpar = x[iParam+1] ;
903 Double_t epar = x[iParam+2] ;
9a2cdbdf 904// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
88cb7938 905 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
a978bcea 906// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
907 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
908//DP: No incident angle dependence in gradient yet!!!!!!
88cb7938 909 Double_t r4 = dr*dr*dr*dr ;
910 Double_t r295 = TMath::Power(dr,2.95) ;
911 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
912 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
913
914 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
915 iParam++ ;
916 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
917 iParam++ ;
918 Grad[iParam] += shape ; // Derivative over energy
919 iParam++ ;
9a1398dd 920 }
921 }
922 efit = 0;
923 iparam = 0 ;
924
925 while(iparam < nPar ){
926 Double_t xpar = x[iparam] ;
927 Double_t zpar = x[iparam+1] ;
928 Double_t epar = x[iparam+2] ;
929 iparam += 3 ;
9a2cdbdf 930// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
a978bcea 931// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
932 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 933 }
934
88cb7938 935 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
9a1398dd 936 // Here we assume, that sigma = sqrt(E)
937 }
83974468 938
d15a28e7 939}
940
941//____________________________________________________________________________
702ab87e 942void AliPHOSClusterizerv1::Print(const Option_t *)const
d15a28e7 943{
d1de15f5 944 // Print clusterizer parameters
945
21cd0c07 946 TString message ;
947 TString taskName(GetName()) ;
948 taskName.ReplaceAll(Version(), "") ;
949
950 if( strcmp(GetName(), "") !=0 ) {
9a1398dd 951 // Print parameters
21cd0c07 952 message = "\n--------------- %s %s -----------\n" ;
953 message += "Clusterizing digits from the file: %s\n" ;
954 message += " Branch: %s\n" ;
955 message += " EMC Clustering threshold = %f\n" ;
956 message += " EMC Local Maximum cut = %f\n" ;
957 message += " EMC Logarothmic weight = %f\n" ;
958 message += " CPV Clustering threshold = %f\n" ;
959 message += " CPV Local Maximum cut = %f\n" ;
960 message += " CPV Logarothmic weight = %f\n" ;
9a1398dd 961 if(fToUnfold)
21cd0c07 962 message += " Unfolding on\n" ;
9a1398dd 963 else
21cd0c07 964 message += " Unfolding off\n" ;
9a1398dd 965
21cd0c07 966 message += "------------------------------------------------------------------" ;
9a1398dd 967 }
21cd0c07 968 else
969 message = " AliPHOSClusterizerv1 not initialized " ;
970
351dd634 971 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
21cd0c07 972 taskName.Data(),
973 GetTitle(),
974 taskName.Data(),
975 GetName(),
976 fEmcClusteringThreshold,
977 fEmcLocMaxCut,
978 fW0,
979 fCpvClusteringThreshold,
980 fCpvLocMaxCut,
351dd634 981 fW0CPV )) ;
9a1398dd 982}
a978bcea 983//____________________________________________________________________________
a4e98857 984void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
985{
986 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
9a1398dd 987
9a2cdbdf 988 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
989 fEMCRecPoints->GetEntriesFast(),
990 fCPVRecPoints->GetEntriesFast() )) ;
88cb7938 991
9a1398dd 992 if(strstr(option,"all")) {
709e117a 993 printf("\n EMC clusters \n") ;
994 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
9a1398dd 995 Int_t index ;
9a2cdbdf 996 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
997 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
9a1398dd 998 TVector3 locpos;
999 rp->GetLocalPosition(locpos);
9a1398dd 1000 Float_t lambda[2];
1001 rp->GetElipsAxis(lambda);
9a1398dd 1002 Int_t * primaries;
1003 Int_t nprimaries;
1004 primaries = rp->GetPrimaries(nprimaries);
709e117a 1005 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
11f9c5ff 1006 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1007 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
3bf72d32 1008
21cd0c07 1009 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
709e117a 1010 printf("%d ", primaries[iprimary] ) ;
21cd0c07 1011 }
709e117a 1012 printf("\n") ;
1013 }
1014
1015 //Now plot CPV recPoints
1016 printf("\n CPV clusters \n") ;
1017 printf("Index Ene(MeV) Module X Y Z \n") ;
9a2cdbdf 1018 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1019 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
2bb500e5 1020
1021 TVector3 locpos;
1022 rp->GetLocalPosition(locpos);
1023
1024 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1025 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1026 locpos.X(), locpos.Y(), locpos.Z()) ;
88cb7938 1027 }
9a1398dd 1028 }
d15a28e7 1029}
9a1398dd 1030
af0b0570 1031
1032//____________________________________________________________________________
1033void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1034{
1035 //For each EMC rec. point set the distance to the nearest bad crystal.
1036 //Author: Boris Polichtchouk
1037
e68222ce 1038 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1039 AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
af0b0570 1040
af0b0570 1041 Int_t badIds[8000];
e68222ce 1042 fgCalibData->EmcBadChannelIds(badIds);
af0b0570 1043
1044 AliPHOSEmcRecPoint* rp;
1045
1046 TMatrixF gmat;
1047 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1048 TVector3 gposBadChannel; // global position of bad crystal
1049 TVector3 dR;
1050
1051 Float_t dist,minDist;
e250de8f 1052 Int_t relid[4] ;
af0b0570 1053
9a2cdbdf 1054 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1055 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
af0b0570 1056 minDist = 1.e+07;
e68222ce 1057 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
e250de8f 1058 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1059 if(relid[0]!=rp->GetPHOSMod())
1060 continue ;
af0b0570 1061 rp->GetGlobalPosition(gposRecPoint,gmat);
9a2cdbdf 1062 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
af0b0570 1063 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1064 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1065 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1066 dR = gposBadChannel-gposRecPoint;
1067 dist = dR.Mag();
1068 if(dist<minDist) minDist = dist;
1069 }
1070
1071 rp->SetDistanceToBadCrystal(minDist);
1072 }
1073
1074}