Loaders removed from the reconstruction code (C.Cheshkov)
[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 *
20 * $Log$
9a2cdbdf 21 * Revision 1.112 2007/08/22 09:20:50 hristov
22 * Updated QA classes (Yves)
23 *
a5fa6165 24 * Revision 1.111 2007/08/08 12:11:28 kharlov
25 * Protection against uninitialized fQADM
26 *
7ae07b77 27 * Revision 1.110 2007/08/07 14:16:00 kharlov
28 * Quality assurance added (Yves Schutz)
29 *
ddd1a39c 30 * Revision 1.109 2007/07/24 17:20:35 policheh
31 * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
32 * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
33 *
3799bcb5 34 * Revision 1.108 2007/06/18 07:00:51 kharlov
35 * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
36 *
cd79ec76 37 * Revision 1.107 2007/05/25 14:12:26 policheh
38 * Local to tracking CS transformation added for CPV rec. points
39 *
76e42795 40 * Revision 1.106 2007/05/24 13:01:22 policheh
41 * Local to tracking CS transformation invoked for each EMC rec.point
42 *
788ab371 43 * Revision 1.105 2007/05/02 13:41:22 kharlov
44 * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
45 *
f36dfa9a 46 * Revision 1.104 2007/04/27 16:55:53 kharlov
47 * Calibration stops if PHOS CDB objects do not exist
48 *
76d78b7a 49 * Revision 1.103 2007/04/11 11:55:45 policheh
50 * SetDistancesToBadChannels() added.
51 *
af0b0570 52 * Revision 1.102 2007/03/28 19:18:15 kharlov
53 * RecPoints recalculation in TSM removed
54 *
8c1fb709 55 * Revision 1.101 2007/03/06 06:51:27 kharlov
56 * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
57 *
a978bcea 58 * Revision 1.100 2007/01/10 11:07:26 kharlov
59 * Raw digits writing to file (B.Polichtchouk)
60 *
cbddd97f 61 * Revision 1.99 2006/11/07 16:49:51 kharlov
62 * Corrections for next event switch in case of raw data (B.Polichtchouk)
63 *
e1c29d81 64 * Revision 1.98 2006/10/27 17:14:27 kharlov
65 * Introduce AliDebug and AliLog (B.Polichtchouk)
66 *
e1afea68 67 * Revision 1.97 2006/08/29 11:41:19 kharlov
68 * Missing implementation of ctors and = operator are added
69 *
f1133801 70 * Revision 1.96 2006/08/25 16:56:30 kharlov
71 * Compliance with Effective C++
72 *
0378398c 73 * Revision 1.95 2006/08/11 12:36:26 cvetan
74 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
75 *
52c5f046 76 * Revision 1.94 2006/08/07 12:27:49 hristov
77 * Removing obsolete code which affected the event numbering scheme
78 *
4d4750bd 79 * Revision 1.93 2006/08/01 12:20:17 cvetan
80 * 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
81 *
f5eaa851 82 * Revision 1.92 2006/04/29 20:26:46 hristov
83 * Separate EMC and CPV calibration (Yu.Kharlov)
84 *
e95226ae 85 * Revision 1.91 2006/04/22 10:30:17 hristov
86 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
87 *
27a73a5d 88 * Revision 1.90 2006/04/11 15:22:59 hristov
89 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
90 *
28871337 91 * Revision 1.89 2006/03/13 14:05:42 kharlov
92 * Calibration objects for EMC and CPV
93 *
fc6706cb 94 * Revision 1.88 2006/01/11 08:54:52 hristov
95 * Additional protection in case no calibration entry was found
96 *
e521c5a0 97 * Revision 1.87 2005/11/22 08:46:43 kharlov
98 * Updated to new CDB (Boris Polichtchouk)
99 *
ef629623 100 * Revision 1.86 2005/11/14 21:52:43 hristov
101 * Coding conventions
102 *
fdf65bb5 103 * Revision 1.85 2005/09/27 16:08:08 hristov
104 * New version of CDB storage framework (A.Colla)
105 *
9e1ceb13 106 * Revision 1.84 2005/09/21 10:02:47 kharlov
107 * Reading calibration from CDB (Boris Polichtchouk)
108 *
8d6d4fb5 109 * Revision 1.82 2005/09/02 15:43:13 kharlov
110 * Add comments in GetCalibrationParameters and Calibrate
111 *
5c149aeb 112 * Revision 1.81 2005/09/02 14:32:07 kharlov
113 * Calibration of raw data
114 *
44ae287e 115 * Revision 1.80 2005/08/24 15:31:36 kharlov
116 * Setting raw digits flag
117 *
9735c615 118 * Revision 1.79 2005/07/25 15:53:53 kharlov
119 * Read raw data
120 *
3cf4f75f 121 * Revision 1.78 2005/05/28 14:19:04 schutz
122 * Compilation warnings fixed by T.P.
123 *
702ab87e 124 */
125
9a1398dd 126//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
d15a28e7 127//////////////////////////////////////////////////////////////////////////////
9a1398dd 128// Clusterization class. Performs clusterization (collects neighbouring active cells) and
a4e98857 129// unfolds the clusters having several local maxima.
130// Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
9a1398dd 131// PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
f035f6ce 132// parameters including input digits branch title, thresholds etc.)
a4e98857 133// This TTask is normally called from Reconstructioner, but can as well be used in
134// standalone mode.
135// Use Case:
9a2cdbdf 136// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
137// root [1] cl->Digits2Clusters(digitsTree,clusterTree)
138// //finds RecPoints in the current event
a4e98857 139// root [2] cl->SetDigitsBranch("digits2")
f035f6ce 140// //sets another title for Digitis (input) branch
a4e98857 141// root [3] cl->SetRecPointsBranch("recp2")
f035f6ce 142// //sets another title four output branches
a4e98857 143// root [4] cl->SetEmcLocalMaxCut(0.03)
9a1398dd 144// //set clusterization parameters
ed4205d8 145
d15a28e7 146// --- ROOT system ---
147
148#include "TMath.h"
9a1398dd 149#include "TMinuit.h"
150#include "TTree.h"
9a1398dd 151#include "TBenchmark.h"
9a2cdbdf 152#include "TClonesArray.h"
d15a28e7 153
154// --- Standard library ---
155
d15a28e7 156// --- AliRoot header files ---
351dd634 157#include "AliLog.h"
9a2cdbdf 158#include "AliConfig.h"
e957fea8 159#include "AliPHOSGeometry.h"
d15a28e7 160#include "AliPHOSClusterizerv1.h"
88cb7938 161#include "AliPHOSEmcRecPoint.h"
9a1398dd 162#include "AliPHOSCpvRecPoint.h"
d15a28e7 163#include "AliPHOSDigit.h"
9a1398dd 164#include "AliPHOSDigitizer.h"
ba898748 165#include "AliPHOSCalibrationDB.h"
9e1ceb13 166#include "AliCDBManager.h"
eaa5012b 167#include "AliCDBStorage.h"
fdf65bb5 168#include "AliCDBEntry.h"
3799bcb5 169#include "AliPHOSRecoParam.h"
ddd1a39c 170#include "AliPHOSQualAssDataMaker.h"
9a2cdbdf 171#include "AliPHOSCalibData.h"
172#include "AliPHOSReconstructor.h"
d15a28e7 173
174ClassImp(AliPHOSClusterizerv1)
f035f6ce 175
d15a28e7 176//____________________________________________________________________________
0378398c 177AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
178 AliPHOSClusterizer(),
179 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
180 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
181 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
182 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
183 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
184 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
9a2cdbdf 185 fW0CPV(0), fEmcTimeGate(0),
186 fIsOldRCUFormat(0)
d15a28e7 187{
f035f6ce 188 // default ctor (to be used mainly by Streamer)
f035f6ce 189
8d0f3f77 190 InitParameters() ;
92f521a9 191 fDefaultInit = kTRUE ;
9a1398dd 192}
7b7c1533 193
9a1398dd 194//____________________________________________________________________________
9a2cdbdf 195AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
196 AliPHOSClusterizer(geom),
0378398c 197 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
198 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
199 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
200 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
201 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
202 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
9a2cdbdf 203 fW0CPV(0), fEmcTimeGate(0),
204 fIsOldRCUFormat(0)
9a1398dd 205{
a4e98857 206 // ctor with the indication of the file where header Tree and digits Tree are stored
9a1398dd 207
8d0f3f77 208 InitParameters() ;
2bd5457f 209 Init() ;
92f521a9 210 fDefaultInit = kFALSE ;
9a1398dd 211}
fc12304f 212
9688c1dd 213//____________________________________________________________________________
214 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
215{
0bc3b8ed 216 // dtor
217
9688c1dd 218}
9a2cdbdf 219
9a1398dd 220//____________________________________________________________________________
e95226ae 221Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
88cb7938 222{
e95226ae 223 // Convert EMC measured amplitude into real energy.
5c149aeb 224 // Calibration parameters are taken from calibration data base for raw data,
225 // or from digitizer parameters for simulated data.
226
44ae287e 227 if(fCalibData){
228 Int_t relId[4];
9a2cdbdf 229 fGeom->AbsToRelNumbering(absId,relId) ;
44ae287e 230 Int_t module = relId[0];
8d6d4fb5 231 Int_t column = relId[3];
eaa5012b 232 Int_t row = relId[2];
e95226ae 233 if(absId <= fEmcCrystals) { // this is EMC
44ae287e 234 fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
e95226ae 235 return amp*fADCchanelEmc ;
44ae287e 236 }
e95226ae 237 }
238 else{ //simulation
239 if(absId <= fEmcCrystals) // this is EMC
240 return fADCpedestalEmc + amp*fADCchanelEmc ;
241 }
242 return 0;
243}
244
245//____________________________________________________________________________
246Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
247{
248 // Convert digitized CPV amplitude into charge.
249 // Calibration parameters are taken from calibration data base for raw data,
250 // or from digitizer parameters for simulated data.
251
252 if(fCalibData){
253 Int_t relId[4];
9a2cdbdf 254 fGeom->AbsToRelNumbering(absId,relId) ;
e95226ae 255 Int_t module = relId[0];
256 Int_t column = relId[3];
257 Int_t row = relId[2];
258 if(absId > fEmcCrystals) { // this is CPV
fc6706cb 259 fADCchanelCpv = fCalibData->GetADCchannelCpv (module,column,row);
260 fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
261 return fADCpedestalCpv + amp*fADCchanelCpv ;
262 }
ba898748 263 }
264 else{ //simulation
e95226ae 265 if(absId > fEmcCrystals) // this is CPV
ba898748 266 return fADCpedestalCpv+ amp*fADCchanelCpv ;
267 }
e95226ae 268 return 0;
3758d9fc 269}
548f0134 270
3758d9fc 271//____________________________________________________________________________
9a2cdbdf 272void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
a4e98857 273{
9a2cdbdf 274 // Steering method to perform clusterization for one event
275 // The input is the tree with digits.
276 // The output is the tree with clusters.
9a1398dd 277
278 if(strstr(option,"tim"))
7d493c2c 279 gBenchmark->Start("PHOSClusterizer");
9a1398dd 280
eabde521 281 if(strstr(option,"print")) {
88cb7938 282 Print() ;
eabde521 283 return ;
284 }
bed9e3fb 285
eaa5012b 286 GetCalibrationParameters() ;
287
9a2cdbdf 288 MakeClusters() ;
9a1398dd 289
9a2cdbdf 290 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
291 fEMCRecPoints->GetEntries()));
292 if(AliLog::GetGlobalDebugLevel()>1)
293 fEMCRecPoints->Print();
88cb7938 294
9a2cdbdf 295 if(fToUnfold)
296 MakeUnfolding();
ddd1a39c 297
298 //makes the quality assurance data
9a2cdbdf 299 if (GetQualAssDataMaker()) {
300 GetQualAssDataMaker()->SetData(fEMCRecPoints) ;
301 GetQualAssDataMaker()->Exec(AliQualAss::kRECPOINTS) ;
302 GetQualAssDataMaker()->SetData(fCPVRecPoints) ;
303 GetQualAssDataMaker()->Exec(AliQualAss::kRECPOINTS) ;
304 }
7ae07b77 305
9a2cdbdf 306 WriteRecPoints();
7b7c1533 307
9a2cdbdf 308 if(strstr(option,"deb"))
309 PrintRecPoints(option) ;
94de8339 310
9a2cdbdf 311 // PLEASE FIX BY MOVING IT TO ALIRECONSTRUCTION !!!
ddd1a39c 312 //Write the quality assurance data only after the last event
9a2cdbdf 313 // if (GetQualAssDataMaker() && fEventCounter == gime->MaxEvent())
314 // GetQualAssDataMaker()->Finish(AliQualAss::kRECPOINTS) ;
ddd1a39c 315
9a1398dd 316 if(strstr(option,"tim")){
317 gBenchmark->Stop("PHOSClusterizer");
9a2cdbdf 318 AliInfo(Form("took %f seconds for Clusterizing\n",
319 gBenchmark->GetCpuTime("PHOSClusterizer")));
88cb7938 320 }
9a1398dd 321}
322
323//____________________________________________________________________________
a0636361 324Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
d1de15f5 325 Int_t nPar, Float_t * fitparameters) const
9a1398dd 326{
327 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
a4e98857 328 // The initial values for fitting procedure are set equal to the positions of local maxima.
f035f6ce 329 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
9a1398dd 330
88cb7938 331
9a1398dd 332 gMinuit->mncler(); // Reset Minuit's list of paramters
333 gMinuit->SetPrintLevel(-1) ; // No Printout
334 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
335 // To set the address of the minimization function
336
337 TList * toMinuit = new TList();
338 toMinuit->AddAt(emcRP,0) ;
9a2cdbdf 339 toMinuit->AddAt(fDigitsArr,1) ;
340 toMinuit->AddAt(fGeom,2) ;
9a1398dd 341
342 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
343
344 // filling initial values for fit parameters
345 AliPHOSDigit * digit ;
346
347 Int_t ierflg = 0;
348 Int_t index = 0 ;
349 Int_t nDigits = (Int_t) nPar / 3 ;
350
351 Int_t iDigit ;
352
9a1398dd 353 for(iDigit = 0; iDigit < nDigits; iDigit++){
a0636361 354 digit = maxAt[iDigit];
9a1398dd 355
356 Int_t relid[4] ;
7b7c1533 357 Float_t x = 0.;
358 Float_t z = 0.;
9a2cdbdf 359 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
360 fGeom->RelPosInModule(relid, x, z) ;
9a1398dd 361
362 Float_t energy = maxAtEnergy[iDigit] ;
363
364 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
365 index++ ;
366 if(ierflg != 0){
21cd0c07 367 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
9a1398dd 368 return kFALSE;
369 }
370 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
371 index++ ;
372 if(ierflg != 0){
21cd0c07 373 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
9a1398dd 374 return kFALSE;
375 }
376 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
377 index++ ;
378 if(ierflg != 0){
21cd0c07 379 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
9a1398dd 380 return kFALSE;
381 }
382 }
383
384 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
385 // depends on it.
386 Double_t p1 = 1.0 ;
387 Double_t p2 = 0.0 ;
388
389 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
390 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
391 gMinuit->SetMaxIterations(5);
392 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
393
394 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
395
396 if(ierflg == 4){ // Minimum not found
21cd0c07 397 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
9a1398dd 398 return kFALSE ;
399 }
400 for(index = 0; index < nPar; index++){
401 Double_t err ;
402 Double_t val ;
403 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
404 fitparameters[index] = val ;
405 }
406
407 delete toMinuit ;
408 return kTRUE;
c3c187e5 409
d15a28e7 410}
411
412//____________________________________________________________________________
3758d9fc 413void AliPHOSClusterizerv1::GetCalibrationParameters()
414{
5c149aeb 415 // Set calibration parameters:
eaa5012b 416 // if calibration database exists, they are read from database,
76d78b7a 417 // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
5c149aeb 418 //
ef629623 419 // It is a user responsilibity to open CDB before reconstruction, for example:
420 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
fbf811ec 421
28871337 422 fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
f36dfa9a 423 if (fCalibData->GetCalibDataEmc() == 0)
424 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
425 if (fCalibData->GetCalibDataCpv() == 0)
426 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
427
3758d9fc 428}
548f0134 429
3758d9fc 430//____________________________________________________________________________
a4e98857 431void AliPHOSClusterizerv1::Init()
432{
2bd5457f 433 // Make all memory allocations which can not be done in default constructor.
434 // Attach the Clusterizer task to the list of PHOS tasks
ba898748 435
9a2cdbdf 436 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
3758d9fc 437
7b7c1533 438 if(!gMinuit)
88cb7938 439 gMinuit = new TMinuit(100);
fbf811ec 440
9a1398dd 441}
7b7c1533 442
9a1398dd 443//____________________________________________________________________________
8d0f3f77 444void AliPHOSClusterizerv1::InitParameters()
445{
446
447 fNumberOfCpvClusters = 0 ;
448 fNumberOfEmcClusters = 0 ;
3799bcb5 449
450 const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
451 if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
452
453 const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv();
454 if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
455
456 fCpvClusteringThreshold = parCpv->GetClusteringThreshold();
457 fEmcClusteringThreshold = parEmc->GetClusteringThreshold();
8d0f3f77 458
3799bcb5 459 fEmcLocMaxCut = parEmc->GetLocalMaxCut();
460 fCpvLocMaxCut = parCpv->GetLocalMaxCut();
ba898748 461
3799bcb5 462 fEmcMinE = parEmc->GetMinE();
463 fCpvMinE = parCpv->GetMinE();
ba898748 464
3799bcb5 465 fW0 = parEmc->GetLogWeight();
466 fW0CPV = parCpv->GetLogWeight();
8d0f3f77 467
468 fEmcTimeGate = 1.e-8 ;
469
470 fToUnfold = kTRUE ;
88cb7938 471
ba898748 472 fWrite = kTRUE ;
473
44ae287e 474 fCalibData = 0 ;
ae911977 475
f5eaa851 476 fIsOldRCUFormat = kFALSE;
8d0f3f77 477}
478
479//____________________________________________________________________________
9a1398dd 480Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
d15a28e7 481{
b2a60966 482 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
483 // = 1 are neighbour
484 // = 2 are not neighbour but do not continue searching
a4e98857 485 // neighbours are defined as digits having at least a common vertex
b2a60966 486 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
487 // which is compared to a digit (d2) not yet in a cluster
488
d15a28e7 489 Int_t rv = 0 ;
490
d15a28e7 491 Int_t relid1[4] ;
9a2cdbdf 492 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
d15a28e7 493
494 Int_t relid2[4] ;
9a2cdbdf 495 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
d15a28e7 496
9688c1dd 497 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
92862013 498 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
499 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
d15a28e7 500
92862013 501 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
9688c1dd 502 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
d15a28e7 503 rv = 1 ;
504 }
505 else {
506 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
88cb7938 507 rv = 2; // Difference in row numbers is too large to look further
d15a28e7 508 }
509
510 }
511 else {
512
9688c1dd 513 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
d15a28e7 514 rv=2 ;
515
516 }
d72dfbc3 517
d15a28e7 518 return rv ;
519}
ba898748 520//____________________________________________________________________________
e95226ae 521void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
522{
523 // Remove digits with amplitudes below threshold
524
ba898748 525 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
526 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
e95226ae 527 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
528 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) )
ba898748 529 digits->RemoveAt(i) ;
530 }
531 digits->Compress() ;
532 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
533 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
534 digit->SetIndexInList(i) ;
535 }
ba898748 536}
9a2cdbdf 537
d15a28e7 538//____________________________________________________________________________
9a1398dd 539Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
d15a28e7 540{
b2a60966 541 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
542
9f616d61 543 Bool_t rv = kFALSE ;
d15a28e7 544
9a2cdbdf 545 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
d15a28e7 546
2b629790 547 if(digit->GetId() <= nEMC ) rv = kTRUE;
ed4205d8 548
549 return rv ;
550}
551
552//____________________________________________________________________________
9a1398dd 553Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
ed4205d8 554{
fad3e5b9 555 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
ed4205d8 556
557 Bool_t rv = kFALSE ;
88cb7938 558
9a2cdbdf 559 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
ed4205d8 560
2b629790 561 if(digit->GetId() > nEMC ) rv = kTRUE;
d15a28e7 562
563 return rv ;
564}
9a1398dd 565
566//____________________________________________________________________________
88cb7938 567void AliPHOSClusterizerv1::WriteRecPoints()
a4e98857 568{
7b7c1533 569
570 // Creates new branches with given title
a4e98857 571 // fills and writes into TreeR.
88cb7938 572
9a1398dd 573 Int_t index ;
ba898748 574 //Evaluate position, dispersion and other RecPoint properties..
9a2cdbdf 575 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
ba898748 576 for(index = 0; index < nEmc; index++){
cd79ec76 577 AliPHOSEmcRecPoint * rp =
9a2cdbdf 578 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
ba898748 579 rp->Purify(fEmcMinE) ;
a978bcea 580 if(rp->GetMultiplicity()==0){
9a2cdbdf 581 fEMCRecPoints->RemoveAt(index) ;
ba898748 582 delete rp ;
cd79ec76 583 continue;
ba898748 584 }
a978bcea 585
cd79ec76 586 // No vertex is available now, calculate corrections in PID
9a2cdbdf 587 rp->EvalAll(fW0,fDigitsArr) ;
cd79ec76 588 TVector3 fakeVtx(0.,0.,0.) ;
9a2cdbdf 589 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
cd79ec76 590 rp->EvalLocal2TrackingCSTransform();
ba898748 591 }
9a2cdbdf 592 fEMCRecPoints->Compress() ;
593 // fEMCRecPoints->Sort() ; //Can not sort until position is calculated!
594 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
595 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
596 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
ba898748 597 }
fbf811ec 598
af0b0570 599 //For each rec.point set the distance to the nearest bad crystal (BVP)
600 SetDistancesToBadChannels();
601
9a1398dd 602 //Now the same for CPV
9a2cdbdf 603 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
604 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
605 rp->EvalAll(fW0CPV,fDigitsArr) ;
76e42795 606 rp->EvalLocal2TrackingCSTransform();
ba898748 607 }
9a2cdbdf 608// fCPVRecPoints->Sort() ;
fbf811ec 609
9a2cdbdf 610 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
611 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
fbf811ec 612
9a2cdbdf 613 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
9a1398dd 614
ba898748 615 if(fWrite){ //We write TreeR
9a2cdbdf 616 fTreeR->Fill();
ba898748 617 }
9a1398dd 618}
619
620//____________________________________________________________________________
621void AliPHOSClusterizerv1::MakeClusters()
622{
623 // Steering method to construct the clusters stored in a list of Reconstructed Points
624 // A cluster is defined as a list of neighbour digits
d15a28e7 625
ba898748 626 //Remove digits below threshold
9a2cdbdf 627 CleanDigits(fDigitsArr) ;
ba898748 628
9a2cdbdf 629 TClonesArray * digitsC = static_cast<TClonesArray*>( fDigitsArr->Clone() ) ;
7b7c1533 630
d15a28e7 631 // Clusterization starts
9a1398dd 632
7b7c1533 633 TIter nextdigit(digitsC) ;
d15a28e7 634 AliPHOSDigit * digit ;
92862013 635 Bool_t notremoved = kTRUE ;
6ad0bfa0 636
88cb7938 637 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
638
639
9a1398dd 640 AliPHOSRecPoint * clu = 0 ;
c161df70 641
d1de15f5 642 TArrayI clusterdigitslist(1500) ;
d15a28e7 643 Int_t index ;
f2bc1b87 644
e95226ae 645 if (( IsInEmc (digit) &&
646 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
647 ( IsInCpv (digit) &&
648 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
d15a28e7 649 Int_t iDigitInCluster = 0 ;
7b7c1533 650
6ad0bfa0 651 if ( IsInEmc(digit) ) {
88cb7938 652 // start a new EMC RecPoint
9a2cdbdf 653 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
654 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
88cb7938 655
9a2cdbdf 656 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
657 clu = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
e95226ae 658 fNumberOfEmcClusters++ ;
659 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
660 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
661 iDigitInCluster++ ;
88cb7938 662 digitsC->Remove(digit) ;
ca77b032 663
d72dfbc3 664 } else {
88cb7938 665
666 // start a new CPV cluster
9a2cdbdf 667 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
668 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
88cb7938 669
9a2cdbdf 670 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
88cb7938 671
9a2cdbdf 672 clu = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
88cb7938 673 fNumberOfCpvClusters++ ;
e95226ae 674 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
88cb7938 675 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
676 iDigitInCluster++ ;
677 digitsC->Remove(digit) ;
d15a28e7 678 nextdigit.Reset() ;
88cb7938 679
680 // Here we remove remaining EMC digits, which cannot make a cluster
681
92862013 682 if( notremoved ) {
88cb7938 683 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
9f616d61 684 if( IsInEmc(digit) )
88cb7938 685 digitsC->Remove(digit) ;
d15a28e7 686 else
88cb7938 687 break ;
688 }
689 notremoved = kFALSE ;
690 }
691
d15a28e7 692 } // else
693
694 nextdigit.Reset() ;
fad3e5b9 695
d15a28e7 696 AliPHOSDigit * digitN ;
697 index = 0 ;
9f616d61 698 while (index < iDigitInCluster){ // scan over digits already in cluster
9a2cdbdf 699 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
88cb7938 700 index++ ;
701 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
702 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
d15a28e7 703 switch (ineb ) {
9f616d61 704 case 0 : // not a neighbour
88cb7938 705 break ;
706 case 1 : // are neighbours
e95226ae 707 if (IsInEmc (digitN))
708 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
709 else
710 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
711
88cb7938 712 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
713 iDigitInCluster++ ;
714 digitsC->Remove(digitN) ;
715 break ;
9f616d61 716 case 2 : // too far from each other
88cb7938 717 goto endofloop;
718 } // switch
719
720 } // while digitN
721
d15a28e7 722 endofloop: ;
88cb7938 723 nextdigit.Reset() ;
724
d15a28e7 725 } // loop over cluster
ad8cfaf4 726
ad8cfaf4 727 } // energy theshold
ad8cfaf4 728
31aa6d6c 729
d15a28e7 730 } // while digit
731
7b7c1533 732 delete digitsC ;
9688c1dd 733
9a1398dd 734}
735
736//____________________________________________________________________________
a4e98857 737void AliPHOSClusterizerv1::MakeUnfolding()
738{
739 // Unfolds clusters using the shape of an ElectroMagnetic shower
9688c1dd 740 // Performs unfolding of all EMC/CPV clusters
9a1398dd 741
a4e98857 742 // Unfold first EMC clusters
9a1398dd 743 if(fNumberOfEmcClusters > 0){
744
9a2cdbdf 745 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 746
747 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
748 Int_t index ;
749 for(index = 0 ; index < numberofNotUnfolded ; index++){
750
9a2cdbdf 751 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
9a1398dd 752 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 753 break ;
9a1398dd 754
755 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 756 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 757 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 758 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
9a1398dd 759
760 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 761 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
9a2cdbdf 762 fEMCRecPoints->Remove(emcRecPoint);
763 fEMCRecPoints->Compress() ;
88cb7938 764 index-- ;
765 fNumberOfEmcClusters -- ;
766 numberofNotUnfolded-- ;
9a1398dd 767 }
c6bf27f2 768 else{
769 emcRecPoint->SetNExMax(1) ; //Only one local maximum
770 }
9a1398dd 771
772 delete[] maxAt ;
773 delete[] maxAtEnergy ;
774 }
775 }
a4e98857 776 // Unfolding of EMC clusters finished
9a1398dd 777
778
a4e98857 779 // Unfold now CPV clusters
9a1398dd 780 if(fNumberOfCpvClusters > 0){
781
9a2cdbdf 782 Int_t nModulesToUnfold = fGeom->GetNModules() ;
9a1398dd 783
784 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
785 Int_t index ;
786 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
787
9a2cdbdf 788 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
9a1398dd 789
790 if(recPoint->GetPHOSMod()> nModulesToUnfold)
88cb7938 791 break ;
9a1398dd 792
88cb7938 793 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
9a1398dd 794
795 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
a0636361 796 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
9a1398dd 797 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
9a2cdbdf 798 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
9a1398dd 799
800 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
88cb7938 801 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
9a2cdbdf 802 fCPVRecPoints->Remove(emcRecPoint);
803 fCPVRecPoints->Compress() ;
88cb7938 804 index-- ;
805 numberofCpvNotUnfolded-- ;
806 fNumberOfCpvClusters-- ;
9a1398dd 807 }
808
809 delete[] maxAt ;
810 delete[] maxAtEnergy ;
811 }
812 }
813 //Unfolding of Cpv clusters finished
814
815}
816
817//____________________________________________________________________________
a978bcea 818Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
9a1398dd 819{
820 // Shape of the shower (see PHOS TDR)
a4e98857 821 // If you change this function, change also the gradient evaluation in ChiSquare()
9a1398dd 822
a978bcea 823 //for the moment we neglect dependence on the incident angle.
824
825 Double_t r2 = x*x + z*z ;
826 Double_t r4 = r2*r2 ;
827 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
9a1398dd 828 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
829 return shape ;
830}
831
832//____________________________________________________________________________
833void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
88cb7938 834 Int_t nMax,
835 AliPHOSDigit ** maxAt,
836 Float_t * maxAtEnergy)
9a1398dd 837{
838 // Performs the unfolding of a cluster with nMax overlapping showers
839
840 Int_t nPar = 3 * nMax ;
841 Float_t * fitparameters = new Float_t[nPar] ;
842
843 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
844 if( !rv ) {
88cb7938 845 // Fit failed, return and remove cluster
c6bf27f2 846 iniEmc->SetNExMax(-1) ;
9a1398dd 847 delete[] fitparameters ;
848 return ;
849 }
850
851 // create ufolded rec points and fill them with new energy lists
852 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
853 // and later correct this number in acordance with actual energy deposition
854
855 Int_t nDigits = iniEmc->GetMultiplicity() ;
856 Float_t * efit = new Float_t[nDigits] ;
a978bcea 857 Float_t xDigit=0.,zDigit=0. ;
7b7c1533 858 Float_t xpar=0.,zpar=0.,epar=0. ;
9a1398dd 859 Int_t relid[4] ;
7b7c1533 860 AliPHOSDigit * digit = 0 ;
9a1398dd 861 Int_t * emcDigits = iniEmc->GetDigitsList() ;
862
a978bcea 863 TVector3 vIncid ;
864
9a1398dd 865 Int_t iparam ;
866 Int_t iDigit ;
867 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 868 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
869 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
870 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
9a1398dd 871 efit[iDigit] = 0;
872
873 iparam = 0 ;
874 while(iparam < nPar ){
875 xpar = fitparameters[iparam] ;
876 zpar = fitparameters[iparam+1] ;
877 epar = fitparameters[iparam+2] ;
878 iparam += 3 ;
9a2cdbdf 879// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
a978bcea 880// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
881 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 882 }
883 }
884
885
886 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
887 // so that energy deposited in each cell is distributed betwin new clusters proportionally
888 // to its contribution to efit
889
890 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
891 Float_t ratio ;
892
893 iparam = 0 ;
894 while(iparam < nPar ){
895 xpar = fitparameters[iparam] ;
896 zpar = fitparameters[iparam+1] ;
897 epar = fitparameters[iparam+2] ;
898 iparam += 3 ;
9a2cdbdf 899// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
9a1398dd 900
7b7c1533 901 AliPHOSEmcRecPoint * emcRP = 0 ;
9a1398dd 902
9a2cdbdf 903 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
9a1398dd 904
9a2cdbdf 905 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
906 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
9a1398dd 907
9a2cdbdf 908 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
909 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
9a1398dd 910 fNumberOfEmcClusters++ ;
c6bf27f2 911 emcRP->SetNExMax((Int_t)nPar/3) ;
9a1398dd 912 }
9a2cdbdf 913 else{//create new entries in fCPVRecPoints
914 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
915 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
9a1398dd 916
9a2cdbdf 917 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
918 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
9a1398dd 919 fNumberOfCpvClusters++ ;
920 }
921
922 Float_t eDigit ;
923 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
9a2cdbdf 924 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
925 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
926 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
a978bcea 927// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
928 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
9a1398dd 929 eDigit = emcEnergies[iDigit] * ratio ;
930 emcRP->AddDigit( *digit, eDigit ) ;
88cb7938 931 }
9a1398dd 932 }
933
934 delete[] fitparameters ;
935 delete[] efit ;
936
937}
938
939//_____________________________________________________________________________
940void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
941{
a4e98857 942 // Calculates the Chi square for the cluster unfolding minimization
9a1398dd 943 // Number of parameters, Gradient, Chi squared, parameters, what to do
944
88cb7938 945 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
9a1398dd 946
88cb7938 947 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
948 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
9a2cdbdf 949 // A bit buggy way to get an access to the geometry
950 // To be revised!
951 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
952
953// TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
9a1398dd 954
88cb7938 955 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
9a1398dd 956
957 Int_t * emcDigits = emcRP->GetDigitsList() ;
958
7b7c1533 959 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
9a1398dd 960
961 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
88cb7938 962
a978bcea 963// TVector3 vInc ;
9a1398dd 964 fret = 0. ;
965 Int_t iparam ;
966
967 if(iflag == 2)
968 for(iparam = 0 ; iparam < nPar ; iparam++)
969 Grad[iparam] = 0 ; // Will evaluate gradient
970
971 Double_t efit ;
972
973 AliPHOSDigit * digit ;
974 Int_t iDigit ;
975
7b7c1533 976 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
9a1398dd 977
88cb7938 978 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
9a1398dd 979
980 Int_t relid[4] ;
981 Float_t xDigit ;
982 Float_t zDigit ;
983
984 geom->AbsToRelNumbering(digit->GetId(), relid) ;
985
986 geom->RelPosInModule(relid, xDigit, zDigit) ;
987
988 if(iflag == 2){ // calculate gradient
989 Int_t iParam = 0 ;
990 efit = 0 ;
991 while(iParam < nPar ){
a978bcea 992 Double_t dx = (xDigit - x[iParam]) ;
88cb7938 993 iParam++ ;
a978bcea 994 Double_t dz = (zDigit - x[iParam]) ;
88cb7938 995 iParam++ ;
9a2cdbdf 996// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
a978bcea 997// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
998 efit += x[iParam] * ShowerShape(dx,dz) ;
88cb7938 999 iParam++ ;
9a1398dd 1000 }
88cb7938 1001 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
9a1398dd 1002 iParam = 0 ;
1003 while(iParam < nPar ){
88cb7938 1004 Double_t xpar = x[iParam] ;
1005 Double_t zpar = x[iParam+1] ;
1006 Double_t epar = x[iParam+2] ;
9a2cdbdf 1007// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
88cb7938 1008 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
a978bcea 1009// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1010 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1011//DP: No incident angle dependence in gradient yet!!!!!!
88cb7938 1012 Double_t r4 = dr*dr*dr*dr ;
1013 Double_t r295 = TMath::Power(dr,2.95) ;
1014 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1015 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1016
1017 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1018 iParam++ ;
1019 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1020 iParam++ ;
1021 Grad[iParam] += shape ; // Derivative over energy
1022 iParam++ ;
9a1398dd 1023 }
1024 }
1025 efit = 0;
1026 iparam = 0 ;
1027
1028 while(iparam < nPar ){
1029 Double_t xpar = x[iparam] ;
1030 Double_t zpar = x[iparam+1] ;
1031 Double_t epar = x[iparam+2] ;
1032 iparam += 3 ;
9a2cdbdf 1033// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
a978bcea 1034// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1035 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
9a1398dd 1036 }
1037
88cb7938 1038 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
9a1398dd 1039 // Here we assume, that sigma = sqrt(E)
1040 }
83974468 1041
d15a28e7 1042}
1043
1044//____________________________________________________________________________
702ab87e 1045void AliPHOSClusterizerv1::Print(const Option_t *)const
d15a28e7 1046{
d1de15f5 1047 // Print clusterizer parameters
1048
21cd0c07 1049 TString message ;
1050 TString taskName(GetName()) ;
1051 taskName.ReplaceAll(Version(), "") ;
1052
1053 if( strcmp(GetName(), "") !=0 ) {
9a1398dd 1054 // Print parameters
21cd0c07 1055 message = "\n--------------- %s %s -----------\n" ;
1056 message += "Clusterizing digits from the file: %s\n" ;
1057 message += " Branch: %s\n" ;
1058 message += " EMC Clustering threshold = %f\n" ;
1059 message += " EMC Local Maximum cut = %f\n" ;
1060 message += " EMC Logarothmic weight = %f\n" ;
1061 message += " CPV Clustering threshold = %f\n" ;
1062 message += " CPV Local Maximum cut = %f\n" ;
1063 message += " CPV Logarothmic weight = %f\n" ;
9a1398dd 1064 if(fToUnfold)
21cd0c07 1065 message += " Unfolding on\n" ;
9a1398dd 1066 else
21cd0c07 1067 message += " Unfolding off\n" ;
9a1398dd 1068
21cd0c07 1069 message += "------------------------------------------------------------------" ;
9a1398dd 1070 }
21cd0c07 1071 else
1072 message = " AliPHOSClusterizerv1 not initialized " ;
1073
351dd634 1074 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
21cd0c07 1075 taskName.Data(),
1076 GetTitle(),
1077 taskName.Data(),
1078 GetName(),
1079 fEmcClusteringThreshold,
1080 fEmcLocMaxCut,
1081 fW0,
1082 fCpvClusteringThreshold,
1083 fCpvLocMaxCut,
351dd634 1084 fW0CPV )) ;
9a1398dd 1085}
a978bcea 1086//____________________________________________________________________________
1087//void AliPHOSClusterizerv1::GetVertex(void)
1088//{ //Extracts vertex posisition
1089//
1090 //ESD
1091//DP - todo if(){
1092//
1093// }
1094
1095// //MC Generator
1096// if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){
1097// Float_t x,y,z ;
1098// gAlice->Generator()->GetOrigin(x,y,z) ;
1099// fVtx.SetXYZ(x,y,z) ;
1100// return ;
1101// }
1102//
1103// //No any source
1104// fVtx[0]=fVtx[1]=fVtx[2]=0. ;
1105//
1106//}
9a1398dd 1107//____________________________________________________________________________
a4e98857 1108void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1109{
1110 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
9a1398dd 1111
9a2cdbdf 1112 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1113 fEMCRecPoints->GetEntriesFast(),
1114 fCPVRecPoints->GetEntriesFast() )) ;
88cb7938 1115
9a1398dd 1116 if(strstr(option,"all")) {
709e117a 1117 printf("\n EMC clusters \n") ;
1118 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
9a1398dd 1119 Int_t index ;
9a2cdbdf 1120 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1121 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
9a1398dd 1122 TVector3 locpos;
1123 rp->GetLocalPosition(locpos);
9a1398dd 1124 Float_t lambda[2];
1125 rp->GetElipsAxis(lambda);
9a1398dd 1126 Int_t * primaries;
1127 Int_t nprimaries;
1128 primaries = rp->GetPrimaries(nprimaries);
709e117a 1129 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
11f9c5ff 1130 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1131 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
3bf72d32 1132
21cd0c07 1133 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
709e117a 1134 printf("%d ", primaries[iprimary] ) ;
21cd0c07 1135 }
709e117a 1136 printf("\n") ;
1137 }
1138
1139 //Now plot CPV recPoints
1140 printf("\n CPV clusters \n") ;
1141 printf("Index Ene(MeV) Module X Y Z \n") ;
9a2cdbdf 1142 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1143 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
2bb500e5 1144
1145 TVector3 locpos;
1146 rp->GetLocalPosition(locpos);
1147
1148 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1149 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1150 locpos.X(), locpos.Y(), locpos.Z()) ;
88cb7938 1151 }
9a1398dd 1152 }
d15a28e7 1153}
9a1398dd 1154
af0b0570 1155
1156//____________________________________________________________________________
1157void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1158{
1159 //For each EMC rec. point set the distance to the nearest bad crystal.
1160 //Author: Boris Polichtchouk
1161
1162 if(!fCalibData->GetNumOfEmcBadChannels()) return;
1163 AliInfo(Form("%d bad channel(s) found.\n",fCalibData->GetNumOfEmcBadChannels()));
1164
af0b0570 1165 Int_t badIds[8000];
1166 fCalibData->EmcBadChannelIds(badIds);
1167
1168 AliPHOSEmcRecPoint* rp;
1169
1170 TMatrixF gmat;
1171 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1172 TVector3 gposBadChannel; // global position of bad crystal
1173 TVector3 dR;
1174
1175 Float_t dist,minDist;
1176
9a2cdbdf 1177 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1178 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
af0b0570 1179 minDist = 1.e+07;
1180
1181 for(Int_t iBad=0; iBad<fCalibData->GetNumOfEmcBadChannels(); iBad++) {
1182 rp->GetGlobalPosition(gposRecPoint,gmat);
9a2cdbdf 1183 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
af0b0570 1184 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1185 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1186 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1187 dR = gposBadChannel-gposRecPoint;
1188 dist = dR.Mag();
1189 if(dist<minDist) minDist = dist;
1190 }
1191
1192 rp->SetDistanceToBadCrystal(minDist);
1193 }
1194
1195}