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