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