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.103 2007/04/11 11:55:45 policheh
22 * SetDistancesToBadChannels() added.
24 * Revision 1.102 2007/03/28 19:18:15 kharlov
25 * RecPoints recalculation in TSM removed
27 * Revision 1.101 2007/03/06 06:51:27 kharlov
28 * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
30 * Revision 1.100 2007/01/10 11:07:26 kharlov
31 * Raw digits writing to file (B.Polichtchouk)
33 * Revision 1.99 2006/11/07 16:49:51 kharlov
34 * Corrections for next event switch in case of raw data (B.Polichtchouk)
36 * Revision 1.98 2006/10/27 17:14:27 kharlov
37 * Introduce AliDebug and AliLog (B.Polichtchouk)
39 * Revision 1.97 2006/08/29 11:41:19 kharlov
40 * Missing implementation of ctors and = operator are added
42 * Revision 1.96 2006/08/25 16:56:30 kharlov
43 * Compliance with Effective C++
45 * Revision 1.95 2006/08/11 12:36:26 cvetan
46 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
48 * Revision 1.94 2006/08/07 12:27:49 hristov
49 * Removing obsolete code which affected the event numbering scheme
51 * Revision 1.93 2006/08/01 12:20:17 cvetan
52 * 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
54 * Revision 1.92 2006/04/29 20:26:46 hristov
55 * Separate EMC and CPV calibration (Yu.Kharlov)
57 * Revision 1.91 2006/04/22 10:30:17 hristov
58 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
60 * Revision 1.90 2006/04/11 15:22:59 hristov
61 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
63 * Revision 1.89 2006/03/13 14:05:42 kharlov
64 * Calibration objects for EMC and CPV
66 * Revision 1.88 2006/01/11 08:54:52 hristov
67 * Additional protection in case no calibration entry was found
69 * Revision 1.87 2005/11/22 08:46:43 kharlov
70 * Updated to new CDB (Boris Polichtchouk)
72 * Revision 1.86 2005/11/14 21:52:43 hristov
75 * Revision 1.85 2005/09/27 16:08:08 hristov
76 * New version of CDB storage framework (A.Colla)
78 * Revision 1.84 2005/09/21 10:02:47 kharlov
79 * Reading calibration from CDB (Boris Polichtchouk)
81 * Revision 1.82 2005/09/02 15:43:13 kharlov
82 * Add comments in GetCalibrationParameters and Calibrate
84 * Revision 1.81 2005/09/02 14:32:07 kharlov
85 * Calibration of raw data
87 * Revision 1.80 2005/08/24 15:31:36 kharlov
88 * Setting raw digits flag
90 * Revision 1.79 2005/07/25 15:53:53 kharlov
93 * Revision 1.78 2005/05/28 14:19:04 schutz
94 * Compilation warnings fixed by T.P.
98 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
99 //////////////////////////////////////////////////////////////////////////////
100 // Clusterization class. Performs clusterization (collects neighbouring active cells) and
101 // unfolds the clusters having several local maxima.
102 // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
103 // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
104 // parameters including input digits branch title, thresholds etc.)
105 // This TTask is normally called from Reconstructioner, but can as well be used in
108 // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
109 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
110 // // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
111 // // and saves recpoints in branch named "recpointsname" (default = "digitsname")
112 // root [1] cl->ExecuteTask()
113 // //finds RecPoints in all events stored in galice.root
114 // root [2] cl->SetDigitsBranch("digits2")
115 // //sets another title for Digitis (input) branch
116 // root [3] cl->SetRecPointsBranch("recp2")
117 // //sets another title four output branches
118 // root [4] cl->SetEmcLocalMaxCut(0.03)
119 // //set clusterization parameters
120 // root [5] cl->ExecuteTask("deb all time")
121 // //once more finds RecPoints options are
122 // // deb - print number of found rec points
123 // // deb all - print number of found RecPoints and some their characteristics
124 // // time - print benchmarking results
126 // --- ROOT system ---
131 #include "TBenchmark.h"
133 // --- Standard library ---
135 // --- AliRoot header files ---
137 #include "AliRunLoader.h"
138 #include "AliGenerator.h"
139 #include "AliPHOSGetter.h"
140 #include "AliPHOSGeometry.h"
141 #include "AliPHOSClusterizerv1.h"
142 #include "AliPHOSEmcRecPoint.h"
143 #include "AliPHOSCpvRecPoint.h"
144 #include "AliPHOSDigit.h"
145 #include "AliPHOSDigitizer.h"
146 #include "AliPHOSCalibrationDB.h"
147 #include "AliCDBManager.h"
148 #include "AliCDBStorage.h"
149 #include "AliCDBEntry.h"
151 ClassImp(AliPHOSClusterizerv1)
153 //____________________________________________________________________________
154 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
155 AliPHOSClusterizer(),
156 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
157 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
158 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
159 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
160 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
161 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
162 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
165 // default ctor (to be used mainly by Streamer)
168 fDefaultInit = kTRUE ;
171 //____________________________________________________________________________
172 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
173 AliPHOSClusterizer(alirunFileName, eventFolderName),
174 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
175 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
176 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
177 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
178 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
179 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
180 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
183 // ctor with the indication of the file where header Tree and digits Tree are stored
187 fDefaultInit = kFALSE ;
190 //____________________________________________________________________________
191 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
192 AliPHOSClusterizer(obj),
193 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
194 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
195 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
196 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
197 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
198 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
199 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
204 //____________________________________________________________________________
205 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
210 //____________________________________________________________________________
211 const TString AliPHOSClusterizerv1::BranchName() const
216 //____________________________________________________________________________
217 Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
219 // Convert EMC measured amplitude into real energy.
220 // Calibration parameters are taken from calibration data base for raw data,
221 // or from digitizer parameters for simulated data.
225 AliPHOSGetter *gime = AliPHOSGetter::Instance();
226 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
227 Int_t module = relId[0];
228 Int_t column = relId[3];
229 Int_t row = relId[2];
230 if(absId <= fEmcCrystals) { // this is EMC
231 fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
232 return amp*fADCchanelEmc ;
236 if(absId <= fEmcCrystals) // this is EMC
237 return fADCpedestalEmc + amp*fADCchanelEmc ;
242 //____________________________________________________________________________
243 Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
245 // Convert digitized CPV amplitude into charge.
246 // Calibration parameters are taken from calibration data base for raw data,
247 // or from digitizer parameters for simulated data.
251 AliPHOSGetter *gime = AliPHOSGetter::Instance();
252 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
253 Int_t module = relId[0];
254 Int_t column = relId[3];
255 Int_t row = relId[2];
256 if(absId > fEmcCrystals) { // this is CPV
257 fADCchanelCpv = fCalibData->GetADCchannelCpv (module,column,row);
258 fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
259 return fADCpedestalCpv + amp*fADCchanelCpv ;
263 if(absId > fEmcCrystals) // this is CPV
264 return fADCpedestalCpv+ amp*fADCchanelCpv ;
269 //____________________________________________________________________________
270 void AliPHOSClusterizerv1::Exec(Option_t *option)
272 // Steering method to perform clusterization for events
273 // in the range from fFirstEvent to fLastEvent.
274 // This range is optionally set by SetEventRange().
275 // if fLastEvent=-1 (by default), then process events until the end.
277 if(strstr(option,"tim"))
278 gBenchmark->Start("PHOSClusterizer");
280 if(strstr(option,"print")) {
285 GetCalibrationParameters() ;
287 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
289 gime->SetRawDigits(kFALSE);
291 gime->SetRawDigits(kTRUE);
293 if (fLastEvent == -1)
294 fLastEvent = gime->MaxEvent() - 1 ;
296 fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
297 Int_t nEvents = fLastEvent - fFirstEvent + 1;
300 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
302 gime->Event(ievent ,"D"); // Read digits from simulated data
304 AliRunLoader * rl = AliRunLoader::GetRunLoader(gime->PhosLoader()->GetTitle());
305 rl->GetEvent(ievent);
306 gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
308 fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
312 AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
313 gime->EmcRecPoints()->GetEntries(),ievent));
314 if(AliLog::GetGlobalDebugLevel()>1)
315 gime->EmcRecPoints()->Print();
322 if(strstr(option,"deb"))
323 PrintRecPoints(option) ;
325 //increment the total number of recpoints per run
326 fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
327 fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
330 if(fWrite) //do not unload in "on flight" mode
333 if(strstr(option,"tim")){
334 gBenchmark->Stop("PHOSClusterizer");
335 AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
336 gBenchmark->GetCpuTime("PHOSClusterizer"),
337 gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
341 //____________________________________________________________________________
342 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
343 Int_t nPar, Float_t * fitparameters) const
345 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
346 // The initial values for fitting procedure are set equal to the positions of local maxima.
347 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
350 AliPHOSGetter * gime = AliPHOSGetter::Instance();
351 TClonesArray * digits = gime->Digits();
353 gMinuit->mncler(); // Reset Minuit's list of paramters
354 gMinuit->SetPrintLevel(-1) ; // No Printout
355 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
356 // To set the address of the minimization function
358 TList * toMinuit = new TList();
359 toMinuit->AddAt(emcRP,0) ;
360 toMinuit->AddAt(digits,1) ;
362 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
364 // filling initial values for fit parameters
365 AliPHOSDigit * digit ;
369 Int_t nDigits = (Int_t) nPar / 3 ;
373 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
375 for(iDigit = 0; iDigit < nDigits; iDigit++){
376 digit = maxAt[iDigit];
381 geom->AbsToRelNumbering(digit->GetId(), relid) ;
382 geom->RelPosInModule(relid, x, z) ;
384 Float_t energy = maxAtEnergy[iDigit] ;
386 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
389 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
392 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
395 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
398 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
401 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
406 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
411 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
412 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
413 gMinuit->SetMaxIterations(5);
414 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
416 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
418 if(ierflg == 4){ // Minimum not found
419 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
422 for(index = 0; index < nPar; index++){
425 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
426 fitparameters[index] = val ;
434 //____________________________________________________________________________
435 void AliPHOSClusterizerv1::GetCalibrationParameters()
437 // Set calibration parameters:
438 // if calibration database exists, they are read from database,
439 // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
441 // It is a user responsilibity to open CDB before reconstruction, for example:
442 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
444 fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
447 //____________________________________________________________________________
448 void AliPHOSClusterizerv1::Init()
450 // Make all memory allocations which can not be done in default constructor.
451 // Attach the Clusterizer task to the list of PHOS tasks
453 AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
455 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
457 AliPHOSGeometry * geom = gime->PHOSGeometry();
459 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
462 gMinuit = new TMinuit(100);
464 if ( !gime->Clusterizer() ) {
465 gime->PostClusterizer(this);
469 //____________________________________________________________________________
470 void AliPHOSClusterizerv1::InitParameters()
473 fNumberOfCpvClusters = 0 ;
474 fNumberOfEmcClusters = 0 ;
476 fCpvClusteringThreshold = 0.0;
477 fEmcClusteringThreshold = 0.2;
479 fEmcLocMaxCut = 0.03 ;
480 fCpvLocMaxCut = 0.03 ;
488 fEmcTimeGate = 1.e-8 ;
492 fRecPointsInRun = 0 ;
498 SetEventRange(0,-1) ;
500 fIsOldRCUFormat = kFALSE;
503 //____________________________________________________________________________
504 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
506 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
508 // = 2 are not neighbour but do not continue searching
509 // neighbours are defined as digits having at least a common vertex
510 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
511 // which is compared to a digit (d2) not yet in a cluster
513 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
518 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
521 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
523 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
524 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
525 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
527 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
528 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
532 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
533 rv = 2; // Difference in row numbers is too large to look further
539 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
546 //____________________________________________________________________________
547 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
549 // Remove digits with amplitudes below threshold
551 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
552 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
553 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
554 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) )
555 digits->RemoveAt(i) ;
558 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
559 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
560 digit->SetIndexInList(i) ;
563 //Overwrite digits tree
564 AliPHOSGetter* gime = AliPHOSGetter::Instance();
565 TTree * treeD = gime->TreeD();
566 treeD->Branch("PHOS", &digits);
568 gime->WriteDigits("OVERWRITE");
569 gime->PhosLoader()->UnloadDigits() ;
571 //____________________________________________________________________________
572 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
574 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
577 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
579 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
581 if(digit->GetId() <= nEMC ) rv = kTRUE;
586 //____________________________________________________________________________
587 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
589 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
593 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
595 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
597 if(digit->GetId() > nEMC ) rv = kTRUE;
602 //____________________________________________________________________________
603 void AliPHOSClusterizerv1::Unload()
605 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
606 gime->PhosLoader()->UnloadDigits() ;
607 gime->PhosLoader()->UnloadRecPoints() ;
610 //____________________________________________________________________________
611 void AliPHOSClusterizerv1::WriteRecPoints()
614 // Creates new branches with given title
615 // fills and writes into TreeR.
617 AliPHOSGetter * gime = AliPHOSGetter::Instance();
619 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
620 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
621 TClonesArray * digits = gime->Digits() ;
624 //Evaluate position, dispersion and other RecPoint properties..
625 Int_t nEmc = emcRecPoints->GetEntriesFast();
626 for(index = 0; index < nEmc; index++){
627 AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
628 rp->Purify(fEmcMinE) ;
629 if(rp->GetMultiplicity()==0){
630 emcRecPoints->RemoveAt(index) ;
634 // No vertex is available now, calculate cirrections in PID
635 rp->EvalAll(fW0,digits) ;
636 TVector3 fakeVtx(0.,0.,0.) ;
637 rp->EvalAll(fW0,fakeVtx,digits) ;
639 emcRecPoints->Compress() ;
640 // emcRecPoints->Sort() ; //Can not sort until position is calculated!
641 // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
642 for(index = 0; index < emcRecPoints->GetEntries(); index++){
643 dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
646 //For each rec.point set the distance to the nearest bad crystal (BVP)
647 SetDistancesToBadChannels();
649 //Now the same for CPV
650 for(index = 0; index < cpvRecPoints->GetEntries(); index++){
651 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
652 rp->EvalAll(fW0CPV,digits) ;
654 // cpvRecPoints->Sort() ;
656 for(index = 0; index < cpvRecPoints->GetEntries(); index++)
657 dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
659 cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
661 if(fWrite){ //We write TreeR
662 TTree * treeR = gime->TreeR();
664 Int_t bufferSize = 32000 ;
665 Int_t splitlevel = 0 ;
668 TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
669 emcBranch->SetTitle(BranchName());
672 TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
673 cpvBranch->SetTitle(BranchName());
678 gime->WriteRecPoints("OVERWRITE");
679 gime->WriteClusterizer("OVERWRITE");
683 //____________________________________________________________________________
684 void AliPHOSClusterizerv1::MakeClusters()
686 // Steering method to construct the clusters stored in a list of Reconstructed Points
687 // A cluster is defined as a list of neighbour digits
690 AliPHOSGetter * gime = AliPHOSGetter::Instance();
692 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
693 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
695 emcRecPoints->Delete() ;
696 cpvRecPoints->Delete() ;
698 TClonesArray * digits = gime->Digits() ;
700 //Remove digits below threshold
701 CleanDigits(digits) ;
703 TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
705 // Clusterization starts
707 TIter nextdigit(digitsC) ;
708 AliPHOSDigit * digit ;
709 Bool_t notremoved = kTRUE ;
711 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
714 AliPHOSRecPoint * clu = 0 ;
716 TArrayI clusterdigitslist(1500) ;
719 if (( IsInEmc (digit) &&
720 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
722 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
723 Int_t iDigitInCluster = 0 ;
725 if ( IsInEmc(digit) ) {
726 // start a new EMC RecPoint
727 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
728 emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
730 emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
731 clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
732 fNumberOfEmcClusters++ ;
733 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
734 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
736 digitsC->Remove(digit) ;
740 // start a new CPV cluster
741 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
742 cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
744 cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
746 clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
747 fNumberOfCpvClusters++ ;
748 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
749 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
751 digitsC->Remove(digit) ;
754 // Here we remove remaining EMC digits, which cannot make a cluster
757 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
759 digitsC->Remove(digit) ;
763 notremoved = kFALSE ;
770 AliPHOSDigit * digitN ;
772 while (index < iDigitInCluster){ // scan over digits already in cluster
773 digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ;
775 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
776 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
778 case 0 : // not a neighbour
780 case 1 : // are neighbours
781 if (IsInEmc (digitN))
782 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
784 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
786 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
788 digitsC->Remove(digitN) ;
790 case 2 : // too far from each other
799 } // loop over cluster
810 //____________________________________________________________________________
811 void AliPHOSClusterizerv1::MakeUnfolding()
813 // Unfolds clusters using the shape of an ElectroMagnetic shower
814 // Performs unfolding of all EMC/CPV clusters
816 AliPHOSGetter * gime = AliPHOSGetter::Instance();
818 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
820 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
821 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
822 TClonesArray * digits = gime->Digits() ;
824 // Unfold first EMC clusters
825 if(fNumberOfEmcClusters > 0){
827 Int_t nModulesToUnfold = geom->GetNModules() ;
829 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
831 for(index = 0 ; index < numberofNotUnfolded ; index++){
833 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
834 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
837 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
838 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
839 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
840 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
842 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
843 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
844 emcRecPoints->Remove(emcRecPoint);
845 emcRecPoints->Compress() ;
847 fNumberOfEmcClusters -- ;
848 numberofNotUnfolded-- ;
851 emcRecPoint->SetNExMax(1) ; //Only one local maximum
855 delete[] maxAtEnergy ;
858 // Unfolding of EMC clusters finished
861 // Unfold now CPV clusters
862 if(fNumberOfCpvClusters > 0){
864 Int_t nModulesToUnfold = geom->GetNModules() ;
866 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
868 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
870 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
872 if(recPoint->GetPHOSMod()> nModulesToUnfold)
875 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
877 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
878 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
879 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
880 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
882 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
883 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
884 cpvRecPoints->Remove(emcRecPoint);
885 cpvRecPoints->Compress() ;
887 numberofCpvNotUnfolded-- ;
888 fNumberOfCpvClusters-- ;
892 delete[] maxAtEnergy ;
895 //Unfolding of Cpv clusters finished
899 //____________________________________________________________________________
900 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
902 // Shape of the shower (see PHOS TDR)
903 // If you change this function, change also the gradient evaluation in ChiSquare()
905 //for the moment we neglect dependence on the incident angle.
907 Double_t r2 = x*x + z*z ;
908 Double_t r4 = r2*r2 ;
909 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
910 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
914 //____________________________________________________________________________
915 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
917 AliPHOSDigit ** maxAt,
918 Float_t * maxAtEnergy)
920 // Performs the unfolding of a cluster with nMax overlapping showers
922 AliPHOSGetter * gime = AliPHOSGetter::Instance();
924 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
926 const TClonesArray * digits = gime->Digits() ;
927 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
928 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
930 Int_t nPar = 3 * nMax ;
931 Float_t * fitparameters = new Float_t[nPar] ;
933 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
935 // Fit failed, return and remove cluster
936 iniEmc->SetNExMax(-1) ;
937 delete[] fitparameters ;
941 // create ufolded rec points and fill them with new energy lists
942 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
943 // and later correct this number in acordance with actual energy deposition
945 Int_t nDigits = iniEmc->GetMultiplicity() ;
946 Float_t * efit = new Float_t[nDigits] ;
947 Float_t xDigit=0.,zDigit=0. ;
948 Float_t xpar=0.,zpar=0.,epar=0. ;
950 AliPHOSDigit * digit = 0 ;
951 Int_t * emcDigits = iniEmc->GetDigitsList() ;
957 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
958 digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
959 geom->AbsToRelNumbering(digit->GetId(), relid) ;
960 geom->RelPosInModule(relid, xDigit, zDigit) ;
964 while(iparam < nPar ){
965 xpar = fitparameters[iparam] ;
966 zpar = fitparameters[iparam+1] ;
967 epar = fitparameters[iparam+2] ;
969 // geom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
970 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
971 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
976 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
977 // so that energy deposited in each cell is distributed betwin new clusters proportionally
978 // to its contribution to efit
980 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
984 while(iparam < nPar ){
985 xpar = fitparameters[iparam] ;
986 zpar = fitparameters[iparam+1] ;
987 epar = fitparameters[iparam+2] ;
989 // geom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
991 AliPHOSEmcRecPoint * emcRP = 0 ;
993 if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
995 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
996 emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
998 (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
999 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
1000 fNumberOfEmcClusters++ ;
1001 emcRP->SetNExMax((Int_t)nPar/3) ;
1003 else{//create new entries in fCpvRecPoints
1004 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
1005 cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
1007 (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
1008 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
1009 fNumberOfCpvClusters++ ;
1013 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
1014 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
1015 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1016 geom->RelPosInModule(relid, xDigit, zDigit) ;
1017 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
1018 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
1019 eDigit = emcEnergies[iDigit] * ratio ;
1020 emcRP->AddDigit( *digit, eDigit ) ;
1024 delete[] fitparameters ;
1029 //_____________________________________________________________________________
1030 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
1032 // Calculates the Chi square for the cluster unfolding minimization
1033 // Number of parameters, Gradient, Chi squared, parameters, what to do
1035 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
1037 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
1038 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
1039 // TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(2)) ; //Vertex position
1041 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
1043 Int_t * emcDigits = emcRP->GetDigitsList() ;
1045 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
1047 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1049 const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
1055 for(iparam = 0 ; iparam < nPar ; iparam++)
1056 Grad[iparam] = 0 ; // Will evaluate gradient
1060 AliPHOSDigit * digit ;
1063 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1065 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
1071 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1073 geom->RelPosInModule(relid, xDigit, zDigit) ;
1075 if(iflag == 2){ // calculate gradient
1078 while(iParam < nPar ){
1079 Double_t dx = (xDigit - x[iParam]) ;
1081 Double_t dz = (zDigit - x[iParam]) ;
1083 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
1084 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
1085 efit += x[iParam] * ShowerShape(dx,dz) ;
1088 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
1090 while(iParam < nPar ){
1091 Double_t xpar = x[iParam] ;
1092 Double_t zpar = x[iParam+1] ;
1093 Double_t epar = x[iParam+2] ;
1094 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1095 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1096 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1097 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1098 //DP: No incident angle dependence in gradient yet!!!!!!
1099 Double_t r4 = dr*dr*dr*dr ;
1100 Double_t r295 = TMath::Power(dr,2.95) ;
1101 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1102 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1104 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1106 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1108 Grad[iParam] += shape ; // Derivative over energy
1115 while(iparam < nPar ){
1116 Double_t xpar = x[iparam] ;
1117 Double_t zpar = x[iparam+1] ;
1118 Double_t epar = x[iparam+2] ;
1120 // geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1121 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1122 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1125 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1126 // Here we assume, that sigma = sqrt(E)
1131 //____________________________________________________________________________
1132 void AliPHOSClusterizerv1::Print(const Option_t *)const
1134 // Print clusterizer parameters
1137 TString taskName(GetName()) ;
1138 taskName.ReplaceAll(Version(), "") ;
1140 if( strcmp(GetName(), "") !=0 ) {
1142 message = "\n--------------- %s %s -----------\n" ;
1143 message += "Clusterizing digits from the file: %s\n" ;
1144 message += " Branch: %s\n" ;
1145 message += " EMC Clustering threshold = %f\n" ;
1146 message += " EMC Local Maximum cut = %f\n" ;
1147 message += " EMC Logarothmic weight = %f\n" ;
1148 message += " CPV Clustering threshold = %f\n" ;
1149 message += " CPV Local Maximum cut = %f\n" ;
1150 message += " CPV Logarothmic weight = %f\n" ;
1152 message += " Unfolding on\n" ;
1154 message += " Unfolding off\n" ;
1156 message += "------------------------------------------------------------------" ;
1159 message = " AliPHOSClusterizerv1 not initialized " ;
1161 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1166 fEmcClusteringThreshold,
1169 fCpvClusteringThreshold,
1173 //____________________________________________________________________________
1174 //void AliPHOSClusterizerv1::GetVertex(void)
1175 //{ //Extracts vertex posisition
1183 // if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){
1185 // gAlice->Generator()->GetOrigin(x,y,z) ;
1186 // fVtx.SetXYZ(x,y,z) ;
1191 // fVtx[0]=fVtx[1]=fVtx[2]=0. ;
1194 //____________________________________________________________________________
1195 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1197 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1199 AliPHOSGetter * gime = AliPHOSGetter::Instance();
1201 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1202 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1204 AliInfo(Form("\nevent %d \n Found %d EMC RecPoints and %d CPV RecPoints",
1205 gAlice->GetEvNumber(),
1206 emcRecPoints->GetEntriesFast(),
1207 cpvRecPoints->GetEntriesFast() )) ;
1209 fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
1210 fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
1213 if(strstr(option,"all")) {
1214 printf("\n EMC clusters \n") ;
1215 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1217 for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1218 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1220 rp->GetLocalPosition(locpos);
1222 rp->GetElipsAxis(lambda);
1225 primaries = rp->GetPrimaries(nprimaries);
1226 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1227 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1228 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1230 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1231 printf("%d ", primaries[iprimary] ) ;
1236 //Now plot CPV recPoints
1237 printf("\n CPV clusters \n") ;
1238 printf("Index Ene(MeV) Module X Y Z \n") ;
1239 for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1240 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
1243 rp->GetLocalPosition(locpos);
1245 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1246 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1247 locpos.X(), locpos.Y(), locpos.Z()) ;
1253 //____________________________________________________________________________
1254 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1256 //For each EMC rec. point set the distance to the nearest bad crystal.
1257 //Author: Boris Polichtchouk
1259 if(!fCalibData->GetNumOfEmcBadChannels()) return;
1260 AliInfo(Form("%d bad channel(s) found.\n",fCalibData->GetNumOfEmcBadChannels()));
1262 AliPHOSGetter* gime = AliPHOSGetter::Instance();
1263 AliPHOSGeometry* geom = gime->PHOSGeometry();
1265 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1268 fCalibData->EmcBadChannelIds(badIds);
1270 AliPHOSEmcRecPoint* rp;
1273 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1274 TVector3 gposBadChannel; // global position of bad crystal
1277 Float_t dist,minDist;
1279 for(Int_t iRP=0; iRP<emcRecPoints->GetEntries(); iRP++){
1280 rp = (AliPHOSEmcRecPoint*)emcRecPoints->At(iRP);
1283 for(Int_t iBad=0; iBad<fCalibData->GetNumOfEmcBadChannels(); iBad++) {
1284 rp->GetGlobalPosition(gposRecPoint,gmat);
1285 geom->RelPosInAlice(badIds[iBad],gposBadChannel);
1286 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1287 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1288 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1289 dR = gposBadChannel-gposRecPoint;
1291 if(dist<minDist) minDist = dist;
1294 rp->SetDistanceToBadCrystal(minDist);