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