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