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