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