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