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