1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /* History of cvs commits:
21 * Revision 1.102 2007/03/28 19:18:15 kharlov
22 * RecPoints recalculation in TSM removed
24 * Revision 1.101 2007/03/06 06:51:27 kharlov
25 * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
27 * Revision 1.100 2007/01/10 11:07:26 kharlov
28 * Raw digits writing to file (B.Polichtchouk)
30 * Revision 1.99 2006/11/07 16:49:51 kharlov
31 * Corrections for next event switch in case of raw data (B.Polichtchouk)
33 * Revision 1.98 2006/10/27 17:14:27 kharlov
34 * Introduce AliDebug and AliLog (B.Polichtchouk)
36 * Revision 1.97 2006/08/29 11:41:19 kharlov
37 * Missing implementation of ctors and = operator are added
39 * Revision 1.96 2006/08/25 16:56:30 kharlov
40 * Compliance with Effective C++
42 * Revision 1.95 2006/08/11 12:36:26 cvetan
43 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
45 * Revision 1.94 2006/08/07 12:27:49 hristov
46 * Removing obsolete code which affected the event numbering scheme
48 * Revision 1.93 2006/08/01 12:20:17 cvetan
49 * 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
51 * Revision 1.92 2006/04/29 20:26:46 hristov
52 * Separate EMC and CPV calibration (Yu.Kharlov)
54 * Revision 1.91 2006/04/22 10:30:17 hristov
55 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
57 * Revision 1.90 2006/04/11 15:22:59 hristov
58 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
60 * Revision 1.89 2006/03/13 14:05:42 kharlov
61 * Calibration objects for EMC and CPV
63 * Revision 1.88 2006/01/11 08:54:52 hristov
64 * Additional protection in case no calibration entry was found
66 * Revision 1.87 2005/11/22 08:46:43 kharlov
67 * Updated to new CDB (Boris Polichtchouk)
69 * Revision 1.86 2005/11/14 21:52:43 hristov
72 * Revision 1.85 2005/09/27 16:08:08 hristov
73 * New version of CDB storage framework (A.Colla)
75 * Revision 1.84 2005/09/21 10:02:47 kharlov
76 * Reading calibration from CDB (Boris Polichtchouk)
78 * Revision 1.82 2005/09/02 15:43:13 kharlov
79 * Add comments in GetCalibrationParameters and Calibrate
81 * Revision 1.81 2005/09/02 14:32:07 kharlov
82 * Calibration of raw data
84 * Revision 1.80 2005/08/24 15:31:36 kharlov
85 * Setting raw digits flag
87 * Revision 1.79 2005/07/25 15:53:53 kharlov
90 * Revision 1.78 2005/05/28 14:19:04 schutz
91 * Compilation warnings fixed by T.P.
95 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
96 //////////////////////////////////////////////////////////////////////////////
97 // Clusterization class. Performs clusterization (collects neighbouring active cells) and
98 // unfolds the clusters having several local maxima.
99 // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
100 // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
101 // parameters including input digits branch title, thresholds etc.)
102 // This TTask is normally called from Reconstructioner, but can as well be used in
105 // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
106 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
107 // // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
108 // // and saves recpoints in branch named "recpointsname" (default = "digitsname")
109 // root [1] cl->ExecuteTask()
110 // //finds RecPoints in all events stored in galice.root
111 // root [2] cl->SetDigitsBranch("digits2")
112 // //sets another title for Digitis (input) branch
113 // root [3] cl->SetRecPointsBranch("recp2")
114 // //sets another title four output branches
115 // root [4] cl->SetEmcLocalMaxCut(0.03)
116 // //set clusterization parameters
117 // root [5] cl->ExecuteTask("deb all time")
118 // //once more finds RecPoints options are
119 // // deb - print number of found rec points
120 // // deb all - print number of found RecPoints and some their characteristics
121 // // time - print benchmarking results
123 // --- ROOT system ---
128 #include "TBenchmark.h"
130 // --- Standard library ---
132 // --- AliRoot header files ---
134 #include "AliRunLoader.h"
135 #include "AliGenerator.h"
136 #include "AliPHOSGetter.h"
137 #include "AliPHOSGeometry.h"
138 #include "AliPHOSClusterizerv1.h"
139 #include "AliPHOSEmcRecPoint.h"
140 #include "AliPHOSCpvRecPoint.h"
141 #include "AliPHOSDigit.h"
142 #include "AliPHOSDigitizer.h"
143 #include "AliPHOSCalibrationDB.h"
144 #include "AliCDBManager.h"
145 #include "AliCDBStorage.h"
146 #include "AliCDBEntry.h"
148 ClassImp(AliPHOSClusterizerv1)
150 //____________________________________________________________________________
151 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
152 AliPHOSClusterizer(),
153 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
154 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
155 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
156 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
157 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
158 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
159 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
162 // default ctor (to be used mainly by Streamer)
165 fDefaultInit = kTRUE ;
168 //____________________________________________________________________________
169 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
170 AliPHOSClusterizer(alirunFileName, eventFolderName),
171 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
172 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
173 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
174 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
175 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
176 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
177 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
180 // ctor with the indication of the file where header Tree and digits Tree are stored
184 fDefaultInit = kFALSE ;
187 //____________________________________________________________________________
188 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
189 AliPHOSClusterizer(obj),
190 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
191 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
192 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
193 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
194 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
195 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
196 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
201 //____________________________________________________________________________
202 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
207 //____________________________________________________________________________
208 const TString AliPHOSClusterizerv1::BranchName() const
213 //____________________________________________________________________________
214 Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
216 // Convert EMC measured amplitude into real energy.
217 // Calibration parameters are taken from calibration data base for raw data,
218 // or from digitizer parameters for simulated data.
222 AliPHOSGetter *gime = AliPHOSGetter::Instance();
223 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
224 Int_t module = relId[0];
225 Int_t column = relId[3];
226 Int_t row = relId[2];
227 if(absId <= fEmcCrystals) { // this is EMC
228 fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
229 return amp*fADCchanelEmc ;
233 if(absId <= fEmcCrystals) // this is EMC
234 return fADCpedestalEmc + amp*fADCchanelEmc ;
239 //____________________________________________________________________________
240 Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
242 // Convert digitized CPV amplitude into charge.
243 // Calibration parameters are taken from calibration data base for raw data,
244 // or from digitizer parameters for simulated data.
248 AliPHOSGetter *gime = AliPHOSGetter::Instance();
249 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
250 Int_t module = relId[0];
251 Int_t column = relId[3];
252 Int_t row = relId[2];
253 if(absId > fEmcCrystals) { // this is CPV
254 fADCchanelCpv = fCalibData->GetADCchannelCpv (module,column,row);
255 fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
256 return fADCpedestalCpv + amp*fADCchanelCpv ;
260 if(absId > fEmcCrystals) // this is CPV
261 return fADCpedestalCpv+ amp*fADCchanelCpv ;
266 //____________________________________________________________________________
267 void AliPHOSClusterizerv1::Exec(Option_t *option)
269 // Steering method to perform clusterization for events
270 // in the range from fFirstEvent to fLastEvent.
271 // This range is optionally set by SetEventRange().
272 // if fLastEvent=-1 (by default), then process events until the end.
274 if(strstr(option,"tim"))
275 gBenchmark->Start("PHOSClusterizer");
277 if(strstr(option,"print")) {
282 GetCalibrationParameters() ;
284 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
286 gime->SetRawDigits(kFALSE);
288 gime->SetRawDigits(kTRUE);
290 if (fLastEvent == -1)
291 fLastEvent = gime->MaxEvent() - 1 ;
293 fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
294 Int_t nEvents = fLastEvent - fFirstEvent + 1;
297 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
299 gime->Event(ievent ,"D"); // Read digits from simulated data
301 AliRunLoader * rl = AliRunLoader::GetRunLoader(gime->PhosLoader()->GetTitle());
302 rl->GetEvent(ievent);
303 gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
305 fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
309 AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
310 gime->EmcRecPoints()->GetEntries(),ievent));
311 if(AliLog::GetGlobalDebugLevel()>1)
312 gime->EmcRecPoints()->Print();
319 if(strstr(option,"deb"))
320 PrintRecPoints(option) ;
322 //increment the total number of recpoints per run
323 fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
324 fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
327 if(fWrite) //do not unload in "on flight" mode
330 if(strstr(option,"tim")){
331 gBenchmark->Stop("PHOSClusterizer");
332 AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
333 gBenchmark->GetCpuTime("PHOSClusterizer"),
334 gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
338 //____________________________________________________________________________
339 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
340 Int_t nPar, Float_t * fitparameters) const
342 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
343 // The initial values for fitting procedure are set equal to the positions of local maxima.
344 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
347 AliPHOSGetter * gime = AliPHOSGetter::Instance();
348 TClonesArray * digits = gime->Digits();
350 gMinuit->mncler(); // Reset Minuit's list of paramters
351 gMinuit->SetPrintLevel(-1) ; // No Printout
352 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
353 // To set the address of the minimization function
355 TList * toMinuit = new TList();
356 toMinuit->AddAt(emcRP,0) ;
357 toMinuit->AddAt(digits,1) ;
359 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
361 // filling initial values for fit parameters
362 AliPHOSDigit * digit ;
366 Int_t nDigits = (Int_t) nPar / 3 ;
370 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
372 for(iDigit = 0; iDigit < nDigits; iDigit++){
373 digit = maxAt[iDigit];
378 geom->AbsToRelNumbering(digit->GetId(), relid) ;
379 geom->RelPosInModule(relid, x, z) ;
381 Float_t energy = maxAtEnergy[iDigit] ;
383 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
386 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
389 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
392 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
395 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
398 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
403 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
408 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
409 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
410 gMinuit->SetMaxIterations(5);
411 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
413 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
415 if(ierflg == 4){ // Minimum not found
416 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
419 for(index = 0; index < nPar; index++){
422 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
423 fitparameters[index] = val ;
431 //____________________________________________________________________________
432 void AliPHOSClusterizerv1::GetCalibrationParameters()
434 // Set calibration parameters:
435 // if calibration database exists, they are read from database,
436 // otherwise, they are taken from digitizer.
438 // It is a user responsilibity to open CDB before reconstruction, for example:
439 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
441 AliPHOSGetter * gime = AliPHOSGetter::Instance();
442 // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
444 fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
448 if ( !gime->Digitizer() )
449 gime->LoadDigitizer();
450 AliPHOSDigitizer * dig = gime->Digitizer();
451 fADCchanelEmc = dig->GetEMCchannel() ;
452 fADCpedestalEmc = dig->GetEMCpedestal();
454 fADCchanelCpv = dig->GetCPVchannel() ;
455 fADCpedestalCpv = dig->GetCPVpedestal() ;
459 //____________________________________________________________________________
460 void AliPHOSClusterizerv1::Init()
462 // Make all memory allocations which can not be done in default constructor.
463 // Attach the Clusterizer task to the list of PHOS tasks
465 AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
467 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
469 AliPHOSGeometry * geom = gime->PHOSGeometry();
471 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
474 gMinuit = new TMinuit(100);
476 if ( !gime->Clusterizer() ) {
477 gime->PostClusterizer(this);
481 //____________________________________________________________________________
482 void AliPHOSClusterizerv1::InitParameters()
485 fNumberOfCpvClusters = 0 ;
486 fNumberOfEmcClusters = 0 ;
488 fCpvClusteringThreshold = 0.0;
489 fEmcClusteringThreshold = 0.2;
491 fEmcLocMaxCut = 0.03 ;
492 fCpvLocMaxCut = 0.03 ;
500 fEmcTimeGate = 1.e-8 ;
504 fRecPointsInRun = 0 ;
510 SetEventRange(0,-1) ;
512 fIsOldRCUFormat = kFALSE;
515 //____________________________________________________________________________
516 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
518 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
520 // = 2 are not neighbour but do not continue searching
521 // neighbours are defined as digits having at least a common vertex
522 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
523 // which is compared to a digit (d2) not yet in a cluster
525 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
530 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
533 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
535 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
536 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
537 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
539 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
540 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
544 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
545 rv = 2; // Difference in row numbers is too large to look further
551 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
558 //____________________________________________________________________________
559 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
561 // Remove digits with amplitudes below threshold
563 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
564 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
565 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
566 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) )
567 digits->RemoveAt(i) ;
570 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
571 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
572 digit->SetIndexInList(i) ;
575 //Overwrite digits tree
576 AliPHOSGetter* gime = AliPHOSGetter::Instance();
577 TTree * treeD = gime->TreeD();
578 treeD->Branch("PHOS", &digits);
580 gime->WriteDigits("OVERWRITE");
581 gime->PhosLoader()->UnloadDigits() ;
583 //____________________________________________________________________________
584 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
586 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
589 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
591 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
593 if(digit->GetId() <= nEMC ) rv = kTRUE;
598 //____________________________________________________________________________
599 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
601 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
605 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
607 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
609 if(digit->GetId() > nEMC ) rv = kTRUE;
614 //____________________________________________________________________________
615 void AliPHOSClusterizerv1::Unload()
617 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
618 gime->PhosLoader()->UnloadDigits() ;
619 gime->PhosLoader()->UnloadRecPoints() ;
622 //____________________________________________________________________________
623 void AliPHOSClusterizerv1::WriteRecPoints()
626 // Creates new branches with given title
627 // fills and writes into TreeR.
629 AliPHOSGetter * gime = AliPHOSGetter::Instance();
631 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
632 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
633 TClonesArray * digits = gime->Digits() ;
636 //Evaluate position, dispersion and other RecPoint properties..
637 Int_t nEmc = emcRecPoints->GetEntriesFast();
638 for(index = 0; index < nEmc; index++){
639 AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
640 rp->Purify(fEmcMinE) ;
641 if(rp->GetMultiplicity()==0){
642 emcRecPoints->RemoveAt(index) ;
646 // No vertex is available now, calculate cirrections in PID
647 rp->EvalAll(fW0,digits) ;
648 TVector3 fakeVtx(0.,0.,0.) ;
649 rp->EvalAll(fW0,fakeVtx,digits) ;
651 emcRecPoints->Compress() ;
652 // emcRecPoints->Sort() ; //Can not sort until position is calculated!
653 // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
654 for(index = 0; index < emcRecPoints->GetEntries(); index++){
655 dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
658 //For each rec.point set the distance to the nearest bad crystal (BVP)
659 SetDistancesToBadChannels();
661 //Now the same for CPV
662 for(index = 0; index < cpvRecPoints->GetEntries(); index++){
663 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
664 rp->EvalAll(fW0CPV,digits) ;
666 // cpvRecPoints->Sort() ;
668 for(index = 0; index < cpvRecPoints->GetEntries(); index++)
669 dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
671 cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
673 if(fWrite){ //We write TreeR
674 TTree * treeR = gime->TreeR();
676 Int_t bufferSize = 32000 ;
677 Int_t splitlevel = 0 ;
680 TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
681 emcBranch->SetTitle(BranchName());
684 TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
685 cpvBranch->SetTitle(BranchName());
690 gime->WriteRecPoints("OVERWRITE");
691 gime->WriteClusterizer("OVERWRITE");
695 //____________________________________________________________________________
696 void AliPHOSClusterizerv1::MakeClusters()
698 // Steering method to construct the clusters stored in a list of Reconstructed Points
699 // A cluster is defined as a list of neighbour digits
702 AliPHOSGetter * gime = AliPHOSGetter::Instance();
704 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
705 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
707 emcRecPoints->Delete() ;
708 cpvRecPoints->Delete() ;
710 TClonesArray * digits = gime->Digits() ;
712 //Remove digits below threshold
713 CleanDigits(digits) ;
715 TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
717 // Clusterization starts
719 TIter nextdigit(digitsC) ;
720 AliPHOSDigit * digit ;
721 Bool_t notremoved = kTRUE ;
723 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
726 AliPHOSRecPoint * clu = 0 ;
728 TArrayI clusterdigitslist(1500) ;
731 if (( IsInEmc (digit) &&
732 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
734 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
735 Int_t iDigitInCluster = 0 ;
737 if ( IsInEmc(digit) ) {
738 // start a new EMC RecPoint
739 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
740 emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
742 emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
743 clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
744 fNumberOfEmcClusters++ ;
745 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
746 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
748 digitsC->Remove(digit) ;
752 // start a new CPV cluster
753 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
754 cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
756 cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
758 clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
759 fNumberOfCpvClusters++ ;
760 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
761 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
763 digitsC->Remove(digit) ;
766 // Here we remove remaining EMC digits, which cannot make a cluster
769 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
771 digitsC->Remove(digit) ;
775 notremoved = kFALSE ;
782 AliPHOSDigit * digitN ;
784 while (index < iDigitInCluster){ // scan over digits already in cluster
785 digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ;
787 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
788 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
790 case 0 : // not a neighbour
792 case 1 : // are neighbours
793 if (IsInEmc (digitN))
794 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
796 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
798 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
800 digitsC->Remove(digitN) ;
802 case 2 : // too far from each other
811 } // loop over cluster
822 //____________________________________________________________________________
823 void AliPHOSClusterizerv1::MakeUnfolding()
825 // Unfolds clusters using the shape of an ElectroMagnetic shower
826 // Performs unfolding of all EMC/CPV clusters
828 AliPHOSGetter * gime = AliPHOSGetter::Instance();
830 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
832 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
833 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
834 TClonesArray * digits = gime->Digits() ;
836 // Unfold first EMC clusters
837 if(fNumberOfEmcClusters > 0){
839 Int_t nModulesToUnfold = geom->GetNModules() ;
841 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
843 for(index = 0 ; index < numberofNotUnfolded ; index++){
845 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
846 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
849 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
850 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
851 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
852 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
854 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
855 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
856 emcRecPoints->Remove(emcRecPoint);
857 emcRecPoints->Compress() ;
859 fNumberOfEmcClusters -- ;
860 numberofNotUnfolded-- ;
863 emcRecPoint->SetNExMax(1) ; //Only one local maximum
867 delete[] maxAtEnergy ;
870 // Unfolding of EMC clusters finished
873 // Unfold now CPV clusters
874 if(fNumberOfCpvClusters > 0){
876 Int_t nModulesToUnfold = geom->GetNModules() ;
878 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
880 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
882 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
884 if(recPoint->GetPHOSMod()> nModulesToUnfold)
887 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
889 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
890 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
891 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
892 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
894 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
895 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
896 cpvRecPoints->Remove(emcRecPoint);
897 cpvRecPoints->Compress() ;
899 numberofCpvNotUnfolded-- ;
900 fNumberOfCpvClusters-- ;
904 delete[] maxAtEnergy ;
907 //Unfolding of Cpv clusters finished
911 //____________________________________________________________________________
912 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
914 // Shape of the shower (see PHOS TDR)
915 // If you change this function, change also the gradient evaluation in ChiSquare()
917 //for the moment we neglect dependence on the incident angle.
919 Double_t r2 = x*x + z*z ;
920 Double_t r4 = r2*r2 ;
921 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
922 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
926 //____________________________________________________________________________
927 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
929 AliPHOSDigit ** maxAt,
930 Float_t * maxAtEnergy)
932 // Performs the unfolding of a cluster with nMax overlapping showers
934 AliPHOSGetter * gime = AliPHOSGetter::Instance();
936 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
938 const TClonesArray * digits = gime->Digits() ;
939 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
940 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
942 Int_t nPar = 3 * nMax ;
943 Float_t * fitparameters = new Float_t[nPar] ;
945 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
947 // Fit failed, return and remove cluster
948 iniEmc->SetNExMax(-1) ;
949 delete[] fitparameters ;
953 // create ufolded rec points and fill them with new energy lists
954 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
955 // and later correct this number in acordance with actual energy deposition
957 Int_t nDigits = iniEmc->GetMultiplicity() ;
958 Float_t * efit = new Float_t[nDigits] ;
959 Float_t xDigit=0.,zDigit=0. ;
960 Float_t xpar=0.,zpar=0.,epar=0. ;
962 AliPHOSDigit * digit = 0 ;
963 Int_t * emcDigits = iniEmc->GetDigitsList() ;
969 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
970 digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
971 geom->AbsToRelNumbering(digit->GetId(), relid) ;
972 geom->RelPosInModule(relid, xDigit, zDigit) ;
976 while(iparam < nPar ){
977 xpar = fitparameters[iparam] ;
978 zpar = fitparameters[iparam+1] ;
979 epar = fitparameters[iparam+2] ;
981 // geom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
982 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
983 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
988 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
989 // so that energy deposited in each cell is distributed betwin new clusters proportionally
990 // to its contribution to efit
992 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
996 while(iparam < nPar ){
997 xpar = fitparameters[iparam] ;
998 zpar = fitparameters[iparam+1] ;
999 epar = fitparameters[iparam+2] ;
1001 // geom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
1003 AliPHOSEmcRecPoint * emcRP = 0 ;
1005 if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
1007 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
1008 emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
1010 (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
1011 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
1012 fNumberOfEmcClusters++ ;
1013 emcRP->SetNExMax((Int_t)nPar/3) ;
1015 else{//create new entries in fCpvRecPoints
1016 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
1017 cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
1019 (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
1020 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
1021 fNumberOfCpvClusters++ ;
1025 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
1026 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
1027 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1028 geom->RelPosInModule(relid, xDigit, zDigit) ;
1029 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
1030 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
1031 eDigit = emcEnergies[iDigit] * ratio ;
1032 emcRP->AddDigit( *digit, eDigit ) ;
1036 delete[] fitparameters ;
1041 //_____________________________________________________________________________
1042 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
1044 // Calculates the Chi square for the cluster unfolding minimization
1045 // Number of parameters, Gradient, Chi squared, parameters, what to do
1047 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
1049 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
1050 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
1051 // TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(2)) ; //Vertex position
1053 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
1055 Int_t * emcDigits = emcRP->GetDigitsList() ;
1057 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
1059 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1061 const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
1067 for(iparam = 0 ; iparam < nPar ; iparam++)
1068 Grad[iparam] = 0 ; // Will evaluate gradient
1072 AliPHOSDigit * digit ;
1075 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1077 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
1083 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1085 geom->RelPosInModule(relid, xDigit, zDigit) ;
1087 if(iflag == 2){ // calculate gradient
1090 while(iParam < nPar ){
1091 Double_t dx = (xDigit - x[iParam]) ;
1093 Double_t dz = (zDigit - x[iParam]) ;
1095 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
1096 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
1097 efit += x[iParam] * ShowerShape(dx,dz) ;
1100 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
1102 while(iParam < nPar ){
1103 Double_t xpar = x[iParam] ;
1104 Double_t zpar = x[iParam+1] ;
1105 Double_t epar = x[iParam+2] ;
1106 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1107 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1108 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1109 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1110 //DP: No incident angle dependence in gradient yet!!!!!!
1111 Double_t r4 = dr*dr*dr*dr ;
1112 Double_t r295 = TMath::Power(dr,2.95) ;
1113 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1114 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1116 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1118 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1120 Grad[iParam] += shape ; // Derivative over energy
1127 while(iparam < nPar ){
1128 Double_t xpar = x[iparam] ;
1129 Double_t zpar = x[iparam+1] ;
1130 Double_t epar = x[iparam+2] ;
1132 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1133 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1134 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1137 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1138 // Here we assume, that sigma = sqrt(E)
1143 //____________________________________________________________________________
1144 void AliPHOSClusterizerv1::Print(const Option_t *)const
1146 // Print clusterizer parameters
1149 TString taskName(GetName()) ;
1150 taskName.ReplaceAll(Version(), "") ;
1152 if( strcmp(GetName(), "") !=0 ) {
1154 message = "\n--------------- %s %s -----------\n" ;
1155 message += "Clusterizing digits from the file: %s\n" ;
1156 message += " Branch: %s\n" ;
1157 message += " EMC Clustering threshold = %f\n" ;
1158 message += " EMC Local Maximum cut = %f\n" ;
1159 message += " EMC Logarothmic weight = %f\n" ;
1160 message += " CPV Clustering threshold = %f\n" ;
1161 message += " CPV Local Maximum cut = %f\n" ;
1162 message += " CPV Logarothmic weight = %f\n" ;
1164 message += " Unfolding on\n" ;
1166 message += " Unfolding off\n" ;
1168 message += "------------------------------------------------------------------" ;
1171 message = " AliPHOSClusterizerv1 not initialized " ;
1173 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1178 fEmcClusteringThreshold,
1181 fCpvClusteringThreshold,
1185 //____________________________________________________________________________
1186 //void AliPHOSClusterizerv1::GetVertex(void)
1187 //{ //Extracts vertex posisition
1195 // if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){
1197 // gAlice->Generator()->GetOrigin(x,y,z) ;
1198 // fVtx.SetXYZ(x,y,z) ;
1203 // fVtx[0]=fVtx[1]=fVtx[2]=0. ;
1206 //____________________________________________________________________________
1207 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1209 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1211 AliPHOSGetter * gime = AliPHOSGetter::Instance();
1213 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1214 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1216 AliInfo(Form("\nevent %d \n Found %d EMC RecPoints and %d CPV RecPoints",
1217 gAlice->GetEvNumber(),
1218 emcRecPoints->GetEntriesFast(),
1219 cpvRecPoints->GetEntriesFast() )) ;
1221 fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
1222 fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
1225 if(strstr(option,"all")) {
1226 printf("\n EMC clusters \n") ;
1227 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1229 for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1230 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1232 rp->GetLocalPosition(locpos);
1234 rp->GetElipsAxis(lambda);
1237 primaries = rp->GetPrimaries(nprimaries);
1238 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1239 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1240 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1242 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1243 printf("%d ", primaries[iprimary] ) ;
1248 //Now plot CPV recPoints
1249 printf("\n CPV clusters \n") ;
1250 printf("Index Ene(MeV) Module X Y Z \n") ;
1251 for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1252 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
1255 rp->GetLocalPosition(locpos);
1257 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1258 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1259 locpos.X(), locpos.Y(), locpos.Z()) ;
1265 //____________________________________________________________________________
1266 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1268 //For each EMC rec. point set the distance to the nearest bad crystal.
1269 //Author: Boris Polichtchouk
1271 if(!fCalibData->GetNumOfEmcBadChannels()) return;
1272 AliInfo(Form("%d bad channel(s) found.\n",fCalibData->GetNumOfEmcBadChannels()));
1274 AliPHOSGetter* gime = AliPHOSGetter::Instance();
1275 AliPHOSGeometry* geom = gime->PHOSGeometry();
1277 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1280 fCalibData->EmcBadChannelIds(badIds);
1282 AliPHOSEmcRecPoint* rp;
1285 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1286 TVector3 gposBadChannel; // global position of bad crystal
1289 Float_t dist,minDist;
1291 for(Int_t iRP=0; iRP<emcRecPoints->GetEntries(); iRP++){
1292 rp = (AliPHOSEmcRecPoint*)emcRecPoints->At(iRP);
1295 for(Int_t iBad=0; iBad<fCalibData->GetNumOfEmcBadChannels(); iBad++) {
1296 rp->GetGlobalPosition(gposRecPoint,gmat);
1297 geom->RelPosInAlice(badIds[iBad],gposBadChannel);
1298 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1299 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1300 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1301 dR = gposBadChannel-gposRecPoint;
1303 if(dist<minDist) minDist = dist;
1306 rp->SetDistanceToBadCrystal(minDist);