Coding rule violations fixed.
[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),
d55c9739 201 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
d15a28e7 202{
f035f6ce 203 // default ctor (to be used mainly by Streamer)
f035f6ce 204
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),
d07a44a8 212 fWrite(0),
213 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
a28333c5 214 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
0378398c 215 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
d55c9739 216 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
9a1398dd 217{
a4e98857 218 // ctor with the indication of the file where header Tree and digits Tree are stored
9a1398dd 219
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 }
6483babc 265 fEMCRecPoints->Delete();
266 fCPVRecPoints->Delete();
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
3758d9fc 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
9a1398dd 377//____________________________________________________________________________
8d0f3f77 378void AliPHOSClusterizerv1::InitParameters()
379{
380
381 fNumberOfCpvClusters = 0 ;
382 fNumberOfEmcClusters = 0 ;
3799bcb5 383
7e88424f 384 const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
385 if(!recoParam) AliFatal("Reconstruction parameters are not set!");
3799bcb5 386
0428f331 387 recoParam->Print();
388
389 fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
390 fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
8d0f3f77 391
7e88424f 392 fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
393 fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
ba898748 394
7e88424f 395 fW0 = recoParam->GetEMCLogWeight();
396 fW0CPV = recoParam->GetCPVLogWeight();
8d0f3f77 397
454c62a4 398 fEmcTimeGate = 1.e-6 ;
7e88424f 399 fEcoreRadius = recoParam->GetEMCEcoreRadius();
8d0f3f77 400
7e88424f 401 fToUnfold = recoParam->EMCToUnfold() ;
88cb7938 402
ba898748 403 fWrite = kTRUE ;
8d0f3f77 404}
405
406//____________________________________________________________________________
9a1398dd 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}
ba898748 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
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
479//____________________________________________________________________________
88cb7938 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();
7e88424f 489 Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
e347da7b 490 TVector3 fakeVtx(0.,0.,0.) ;
ba898748 491 for(index = 0; index < nEmc; index++){
cd79ec76 492 AliPHOSEmcRecPoint * rp =
9a2cdbdf 493 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
a28333c5 494 rp->Purify(emcMinE) ;
a978bcea 495 if(rp->GetMultiplicity()==0){
9a2cdbdf 496 fEMCRecPoints->RemoveAt(index) ;
ba898748 497 delete rp ;
cd79ec76 498 continue;
ba898748 499 }
a978bcea 500
cd79ec76 501 // No vertex is available now, calculate corrections in PID
e347da7b 502 rp->EvalAll(fDigitsArr) ;
d55c9739 503 rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
9a2cdbdf 504 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
cd79ec76 505 rp->EvalLocal2TrackingCSTransform();
ba898748 506 }
9a2cdbdf 507 fEMCRecPoints->Compress() ;
6d0a5c2e 508 fEMCRecPoints->Sort() ;
9a2cdbdf 509 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
510 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
511 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
ba898748 512 }
fbf811ec 513
af0b0570 514 //For each rec.point set the distance to the nearest bad crystal (BVP)
515 SetDistancesToBadChannels();
516
9a1398dd 517 //Now the same for CPV
9a2cdbdf 518 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
519 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
e347da7b 520 rp->EvalAll(fDigitsArr) ;
521 rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
76e42795 522 rp->EvalLocal2TrackingCSTransform();
ba898748 523 }
6d0a5c2e 524 fCPVRecPoints->Sort() ;
fbf811ec 525
9a2cdbdf 526 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
527 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
fbf811ec 528
9a2cdbdf 529 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
9a1398dd 530
ba898748 531 if(fWrite){ //We write TreeR
9a2cdbdf 532 fTreeR->Fill();
ba898748 533 }
9a1398dd 534}
535
536//____________________________________________________________________________
537void AliPHOSClusterizerv1::MakeClusters()
538{
539 // Steering method to construct the clusters stored in a list of Reconstructed Points
540 // A cluster is defined as a list of neighbour digits
d15a28e7 541
8d8258f6 542 fNumberOfCpvClusters = 0 ;
543 fNumberOfEmcClusters = 0 ;
544
d07a44a8 545 //Mark all digits as unused yet
78b28884 546 const Int_t maxNDigits = 1500;
d07a44a8 547 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
78b28884 548 if (nDigits > maxNDigits) {
549 AliWarning(Form("Skip event with too high digit occupancy: nDigits=%d",nDigits));
550 return;
551 }
552
d07a44a8 553 for(Int_t i=0; i<nDigits; i++){
554 fDigitsUsed[i]=0 ;
555 }
556 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
557 //e.g. first digit in this module, first CPV digit etc.
d15a28e7 558 AliPHOSDigit * digit ;
78b28884 559 TArrayI clusterdigitslist(maxNDigits) ;
d07a44a8 560 AliPHOSRecPoint * clu = 0 ;
561 for(Int_t i=0; i<nDigits; i++){
562 if(fDigitsUsed[i])
563 continue ;
6ad0bfa0 564
d07a44a8 565 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
88cb7938 566
d07a44a8 567 clu=0 ;
c161df70 568
d15a28e7 569 Int_t index ;
f2bc1b87 570
d07a44a8 571 //is this digit so energetic that start cluster?
a28333c5 572 if (( IsInEmc(digit) && digit->GetEnergy() > fEmcClusteringThreshold ) ||
573 ( IsInCpv(digit) && digit->GetEnergy() > fCpvClusteringThreshold ) ) {
d15a28e7 574 Int_t iDigitInCluster = 0 ;
7b7c1533 575
6ad0bfa0 576 if ( IsInEmc(digit) ) {
88cb7938 577 // start a new EMC RecPoint
9a2cdbdf 578 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
579 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
88cb7938 580
9a2cdbdf 581 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
d07a44a8 582 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
e95226ae 583 fNumberOfEmcClusters++ ;
a28333c5 584 clu->AddDigit(*digit, digit->GetEnergy()) ;
e95226ae 585 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
586 iDigitInCluster++ ;
d07a44a8 587 fDigitsUsed[i]=kTRUE ;
d72dfbc3 588 } else {
88cb7938 589 // start a new CPV cluster
9a2cdbdf 590 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
591 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
88cb7938 592
9a2cdbdf 593 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
d07a44a8 594 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
88cb7938 595 fNumberOfCpvClusters++ ;
a28333c5 596 clu->AddDigit(*digit, digit->GetEnergy()) ;
88cb7938 597 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
598 iDigitInCluster++ ;
d07a44a8 599 fDigitsUsed[i]=kTRUE ;
d15a28e7 600 } // else
601
d07a44a8 602 //Now scan remaining digits in list to find neigbours of our seed
603
d15a28e7 604 AliPHOSDigit * digitN ;
605 index = 0 ;
9f616d61 606 while (index < iDigitInCluster){ // scan over digits already in cluster
d07a44a8 607 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
88cb7938 608 index++ ;
d07a44a8 609 for(Int_t j=iFirst; j<nDigits; j++){
610 if(fDigitsUsed[j])
611 continue ; //look through remaining digits
612 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
88cb7938 613 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
d15a28e7 614 switch (ineb ) {
d07a44a8 615 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
616 iFirst=j ;
617 break ;
9f616d61 618 case 0 : // not a neighbour
88cb7938 619 break ;
620 case 1 : // are neighbours
a28333c5 621 clu->AddDigit(*digitN, digitN->GetEnergy());
d07a44a8 622 clusterdigitslist[iDigitInCluster] = j ;
88cb7938 623 iDigitInCluster++ ;
d07a44a8 624 fDigitsUsed[j]=kTRUE ;
88cb7938 625 break ;
9f616d61 626 case 2 : // too far from each other
88cb7938 627 goto endofloop;
628 } // switch
629
d07a44a8 630 }
88cb7938 631
d07a44a8 632 endofloop: ; //scanned all possible neighbours for this digit
88cb7938 633
d15a28e7 634 } // loop over cluster
ad8cfaf4 635 } // energy theshold
d07a44a8 636 }
9688c1dd 637
9a1398dd 638}
639
640//____________________________________________________________________________
a4e98857 641void AliPHOSClusterizerv1::MakeUnfolding()
642{
643 // Unfolds clusters using the shape of an ElectroMagnetic shower
9688c1dd 644 // Performs unfolding of all EMC/CPV clusters
9a1398dd 645
a4e98857 646 // Unfold first EMC clusters
9a1398dd 647 if(fNumberOfEmcClusters > 0){
648
9a2cdbdf 649 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 650
651 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
652 Int_t index ;
653 for(index = 0 ; index < numberofNotUnfolded ; index++){
654
9a2cdbdf 655 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
9a1398dd 656 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 657 break ;
9a1398dd 658
659 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 660 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 661 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 662 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
9a1398dd 663
664 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 665 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
a28333c5 666
9a2cdbdf 667 fEMCRecPoints->Remove(emcRecPoint);
668 fEMCRecPoints->Compress() ;
88cb7938 669 index-- ;
670 fNumberOfEmcClusters -- ;
671 numberofNotUnfolded-- ;
9a1398dd 672 }
c6bf27f2 673 else{
674 emcRecPoint->SetNExMax(1) ; //Only one local maximum
675 }
9a1398dd 676
677 delete[] maxAt ;
678 delete[] maxAtEnergy ;
679 }
680 }
a4e98857 681 // Unfolding of EMC clusters finished
9a1398dd 682
683
a4e98857 684 // Unfold now CPV clusters
9a1398dd 685 if(fNumberOfCpvClusters > 0){
686
9a2cdbdf 687 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 688
689 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
690 Int_t index ;
691 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
692
9a2cdbdf 693 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
9a1398dd 694
695 if(recPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 696 break ;
9a1398dd 697
88cb7938 698 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
9a1398dd 699
700 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 701 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 702 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 703 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
9a1398dd 704
705 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 706 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
9a2cdbdf 707 fCPVRecPoints->Remove(emcRecPoint);
708 fCPVRecPoints->Compress() ;
88cb7938 709 index-- ;
710 numberofCpvNotUnfolded-- ;
711 fNumberOfCpvClusters-- ;
9a1398dd 712 }
713
714 delete[] maxAt ;
715 delete[] maxAtEnergy ;
716 }
717 }
718 //Unfolding of Cpv clusters finished
719
720}
721
722//____________________________________________________________________________
a978bcea 723Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
9a1398dd 724{
725 // Shape of the shower (see PHOS TDR)
a4e98857 726 // If you change this function, change also the gradient evaluation in ChiSquare()
9a1398dd 727
a978bcea 728 //for the moment we neglect dependence on the incident angle.
729
730 Double_t r2 = x*x + z*z ;
731 Double_t r4 = r2*r2 ;
732 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
9a1398dd 733 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
734 return shape ;
735}
736
737//____________________________________________________________________________
738void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
88cb7938 739 Int_t nMax,
740 AliPHOSDigit ** maxAt,
741 Float_t * maxAtEnergy)
9a1398dd 742{
743 // Performs the unfolding of a cluster with nMax overlapping showers
744
745 Int_t nPar = 3 * nMax ;
746 Float_t * fitparameters = new Float_t[nPar] ;
747
748 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
a28333c5 749
9a1398dd 750 if( !rv ) {
88cb7938 751 // Fit failed, return and remove cluster
c6bf27f2 752 iniEmc->SetNExMax(-1) ;
9a1398dd 753 delete[] fitparameters ;
754 return ;
755 }
756
757 // create ufolded rec points and fill them with new energy lists
758 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
759 // and later correct this number in acordance with actual energy deposition
760
761 Int_t nDigits = iniEmc->GetMultiplicity() ;
762 Float_t * efit = new Float_t[nDigits] ;
a978bcea 763 Float_t xDigit=0.,zDigit=0. ;
7b7c1533 764 Float_t xpar=0.,zpar=0.,epar=0. ;
9a1398dd 765 Int_t relid[4] ;
7b7c1533 766 AliPHOSDigit * digit = 0 ;
9a1398dd 767 Int_t * emcDigits = iniEmc->GetDigitsList() ;
768
a978bcea 769 TVector3 vIncid ;
770
9a1398dd 771 Int_t iparam ;
772 Int_t iDigit ;
773 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 774 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
775 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
776 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
9a1398dd 777 efit[iDigit] = 0;
778
779 iparam = 0 ;
780 while(iparam < nPar ){
781 xpar = fitparameters[iparam] ;
782 zpar = fitparameters[iparam+1] ;
783 epar = fitparameters[iparam+2] ;
784 iparam += 3 ;
9a2cdbdf 785// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
a978bcea 786// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
787 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 788 }
789 }
790
9a1398dd 791 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
792 // so that energy deposited in each cell is distributed betwin new clusters proportionally
793 // to its contribution to efit
794
795 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
796 Float_t ratio ;
797
798 iparam = 0 ;
799 while(iparam < nPar ){
800 xpar = fitparameters[iparam] ;
801 zpar = fitparameters[iparam+1] ;
802 epar = fitparameters[iparam+2] ;
803 iparam += 3 ;
9a2cdbdf 804// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
a28333c5 805
7b7c1533 806 AliPHOSEmcRecPoint * emcRP = 0 ;
9a1398dd 807
9a2cdbdf 808 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
9a1398dd 809
9a2cdbdf 810 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
811 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
9a1398dd 812
9a2cdbdf 813 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
814 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
9a1398dd 815 fNumberOfEmcClusters++ ;
c6bf27f2 816 emcRP->SetNExMax((Int_t)nPar/3) ;
9a1398dd 817 }
9a2cdbdf 818 else{//create new entries in fCPVRecPoints
819 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
820 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
9a1398dd 821
9a2cdbdf 822 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
823 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
9a1398dd 824 fNumberOfCpvClusters++ ;
825 }
826
827 Float_t eDigit ;
828 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 829 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
830 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
831 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
a978bcea 832// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
833 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
9a1398dd 834 eDigit = emcEnergies[iDigit] * ratio ;
835 emcRP->AddDigit( *digit, eDigit ) ;
88cb7938 836 }
9a1398dd 837 }
838
839 delete[] fitparameters ;
840 delete[] efit ;
841
842}
843
844//_____________________________________________________________________________
845void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
846{
a4e98857 847 // Calculates the Chi square for the cluster unfolding minimization
9a1398dd 848 // Number of parameters, Gradient, Chi squared, parameters, what to do
849
88cb7938 850 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
9a1398dd 851
88cb7938 852 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
853 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
9a2cdbdf 854 // A bit buggy way to get an access to the geometry
855 // To be revised!
856 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
857
858// TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
9a1398dd 859
88cb7938 860 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
9a1398dd 861
862 Int_t * emcDigits = emcRP->GetDigitsList() ;
863
7b7c1533 864 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
9a1398dd 865
866 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
88cb7938 867
a978bcea 868// TVector3 vInc ;
9a1398dd 869 fret = 0. ;
870 Int_t iparam ;
871
872 if(iflag == 2)
873 for(iparam = 0 ; iparam < nPar ; iparam++)
874 Grad[iparam] = 0 ; // Will evaluate gradient
875
876 Double_t efit ;
877
878 AliPHOSDigit * digit ;
879 Int_t iDigit ;
880
7b7c1533 881 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
9a1398dd 882
88cb7938 883 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
9a1398dd 884
885 Int_t relid[4] ;
886 Float_t xDigit ;
887 Float_t zDigit ;
888
889 geom->AbsToRelNumbering(digit->GetId(), relid) ;
890
891 geom->RelPosInModule(relid, xDigit, zDigit) ;
892
893 if(iflag == 2){ // calculate gradient
894 Int_t iParam = 0 ;
895 efit = 0 ;
896 while(iParam < nPar ){
a978bcea 897 Double_t dx = (xDigit - x[iParam]) ;
88cb7938 898 iParam++ ;
a978bcea 899 Double_t dz = (zDigit - x[iParam]) ;
88cb7938 900 iParam++ ;
9a2cdbdf 901// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
a978bcea 902// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
903 efit += x[iParam] * ShowerShape(dx,dz) ;
88cb7938 904 iParam++ ;
9a1398dd 905 }
88cb7938 906 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
9a1398dd 907 iParam = 0 ;
908 while(iParam < nPar ){
88cb7938 909 Double_t xpar = x[iParam] ;
910 Double_t zpar = x[iParam+1] ;
911 Double_t epar = x[iParam+2] ;
9a2cdbdf 912// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
88cb7938 913 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
a978bcea 914// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
915 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
916//DP: No incident angle dependence in gradient yet!!!!!!
88cb7938 917 Double_t r4 = dr*dr*dr*dr ;
918 Double_t r295 = TMath::Power(dr,2.95) ;
919 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
920 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
921
922 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
923 iParam++ ;
924 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
925 iParam++ ;
926 Grad[iParam] += shape ; // Derivative over energy
927 iParam++ ;
9a1398dd 928 }
929 }
930 efit = 0;
931 iparam = 0 ;
932
933 while(iparam < nPar ){
934 Double_t xpar = x[iparam] ;
935 Double_t zpar = x[iparam+1] ;
936 Double_t epar = x[iparam+2] ;
937 iparam += 3 ;
9a2cdbdf 938// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
a978bcea 939// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
940 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 941 }
942
88cb7938 943 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
9a1398dd 944 // Here we assume, that sigma = sqrt(E)
945 }
83974468 946
d15a28e7 947}
948
949//____________________________________________________________________________
702ab87e 950void AliPHOSClusterizerv1::Print(const Option_t *)const
d15a28e7 951{
d1de15f5 952 // Print clusterizer parameters
953
21cd0c07 954 TString message ;
955 TString taskName(GetName()) ;
956 taskName.ReplaceAll(Version(), "") ;
957
958 if( strcmp(GetName(), "") !=0 ) {
9a1398dd 959 // Print parameters
21cd0c07 960 message = "\n--------------- %s %s -----------\n" ;
961 message += "Clusterizing digits from the file: %s\n" ;
962 message += " Branch: %s\n" ;
963 message += " EMC Clustering threshold = %f\n" ;
964 message += " EMC Local Maximum cut = %f\n" ;
965 message += " EMC Logarothmic weight = %f\n" ;
966 message += " CPV Clustering threshold = %f\n" ;
967 message += " CPV Local Maximum cut = %f\n" ;
968 message += " CPV Logarothmic weight = %f\n" ;
9a1398dd 969 if(fToUnfold)
21cd0c07 970 message += " Unfolding on\n" ;
9a1398dd 971 else
21cd0c07 972 message += " Unfolding off\n" ;
9a1398dd 973
21cd0c07 974 message += "------------------------------------------------------------------" ;
9a1398dd 975 }
21cd0c07 976 else
977 message = " AliPHOSClusterizerv1 not initialized " ;
978
351dd634 979 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
21cd0c07 980 taskName.Data(),
981 GetTitle(),
982 taskName.Data(),
983 GetName(),
984 fEmcClusteringThreshold,
985 fEmcLocMaxCut,
986 fW0,
987 fCpvClusteringThreshold,
988 fCpvLocMaxCut,
351dd634 989 fW0CPV )) ;
9a1398dd 990}
a978bcea 991//____________________________________________________________________________
a4e98857 992void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
993{
994 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
9a1398dd 995
9a2cdbdf 996 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
997 fEMCRecPoints->GetEntriesFast(),
998 fCPVRecPoints->GetEntriesFast() )) ;
88cb7938 999
9a1398dd 1000 if(strstr(option,"all")) {
709e117a 1001 printf("\n EMC clusters \n") ;
1002 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
9a1398dd 1003 Int_t index ;
9a2cdbdf 1004 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1005 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
9a1398dd 1006 TVector3 locpos;
1007 rp->GetLocalPosition(locpos);
9a1398dd 1008 Float_t lambda[2];
1009 rp->GetElipsAxis(lambda);
9a1398dd 1010 Int_t * primaries;
1011 Int_t nprimaries;
1012 primaries = rp->GetPrimaries(nprimaries);
709e117a 1013 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
11f9c5ff 1014 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1015 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
3bf72d32 1016
21cd0c07 1017 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
709e117a 1018 printf("%d ", primaries[iprimary] ) ;
21cd0c07 1019 }
709e117a 1020 printf("\n") ;
1021 }
1022
1023 //Now plot CPV recPoints
1024 printf("\n CPV clusters \n") ;
1025 printf("Index Ene(MeV) Module X Y Z \n") ;
9a2cdbdf 1026 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1027 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
2bb500e5 1028
1029 TVector3 locpos;
1030 rp->GetLocalPosition(locpos);
1031
1032 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1033 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1034 locpos.X(), locpos.Y(), locpos.Z()) ;
88cb7938 1035 }
9a1398dd 1036 }
d15a28e7 1037}
9a1398dd 1038
af0b0570 1039
1040//____________________________________________________________________________
1041void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1042{
1043 //For each EMC rec. point set the distance to the nearest bad crystal.
1044 //Author: Boris Polichtchouk
1045
e68222ce 1046 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1047 AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
af0b0570 1048
af0b0570 1049 Int_t badIds[8000];
e68222ce 1050 fgCalibData->EmcBadChannelIds(badIds);
af0b0570 1051
1052 AliPHOSEmcRecPoint* rp;
1053
1054 TMatrixF gmat;
1055 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1056 TVector3 gposBadChannel; // global position of bad crystal
1057 TVector3 dR;
1058
1059 Float_t dist,minDist;
79987dcd 1060 Int_t relid[4]={0,0,0,0} ;
1061 TVector3 lpos ;
9a2cdbdf 1062 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1063 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
79987dcd 1064 //evaluate distance to border
1065 relid[0]=rp->GetPHOSMod() ;
1066 relid[2]=1 ;
1067 relid[3]=1 ;
1068 Float_t xcorner,zcorner;
1069 fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1070 rp->GetLocalPosition(lpos) ;
1071 minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
e68222ce 1072 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
e250de8f 1073 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
79987dcd 1074 if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1075 continue ; //bad channels can be in the module which does not exist in simulations.
af0b0570 1076 rp->GetGlobalPosition(gposRecPoint,gmat);
9a2cdbdf 1077 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
af0b0570 1078 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1079 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1080 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1081 dR = gposBadChannel-gposRecPoint;
1082 dist = dR.Mag();
1083 if(dist<minDist) minDist = dist;
1084 }
1085
1086 rp->SetDistanceToBadCrystal(minDist);
1087 }
1088
1089}