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