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