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:
20 * $Log: AliPHOSClusterizerv1.cxx,v $
21 * Revision 1.118 2007/12/11 21:23:26 kharlov
22 * Added possibility to swith off unfolding
24 * Revision 1.117 2007/10/18 08:42:05 kharlov
25 * Bad channels cleaned before clusterization
27 * Revision 1.116 2007/10/01 20:24:08 kharlov
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.
33 * Revision 1.114 2007/09/06 16:06:44 kharlov
34 * Absence of sorting results in loose of all unfolded clusters
36 * Revision 1.113 2007/08/28 12:55:07 policheh
37 * Loaders removed from the reconstruction code (C.Cheshkov)
39 * Revision 1.112 2007/08/22 09:20:50 hristov
40 * Updated QA classes (Yves)
42 * Revision 1.111 2007/08/08 12:11:28 kharlov
43 * Protection against uninitialized fQADM
45 * Revision 1.110 2007/08/07 14:16:00 kharlov
46 * Quality assurance added (Yves Schutz)
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).
52 * Revision 1.108 2007/06/18 07:00:51 kharlov
53 * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
55 * Revision 1.107 2007/05/25 14:12:26 policheh
56 * Local to tracking CS transformation added for CPV rec. points
58 * Revision 1.106 2007/05/24 13:01:22 policheh
59 * Local to tracking CS transformation invoked for each EMC rec.point
61 * Revision 1.105 2007/05/02 13:41:22 kharlov
62 * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
64 * Revision 1.104 2007/04/27 16:55:53 kharlov
65 * Calibration stops if PHOS CDB objects do not exist
67 * Revision 1.103 2007/04/11 11:55:45 policheh
68 * SetDistancesToBadChannels() added.
70 * Revision 1.102 2007/03/28 19:18:15 kharlov
71 * RecPoints recalculation in TSM removed
73 * Revision 1.101 2007/03/06 06:51:27 kharlov
74 * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
76 * Revision 1.100 2007/01/10 11:07:26 kharlov
77 * Raw digits writing to file (B.Polichtchouk)
79 * Revision 1.99 2006/11/07 16:49:51 kharlov
80 * Corrections for next event switch in case of raw data (B.Polichtchouk)
82 * Revision 1.98 2006/10/27 17:14:27 kharlov
83 * Introduce AliDebug and AliLog (B.Polichtchouk)
85 * Revision 1.97 2006/08/29 11:41:19 kharlov
86 * Missing implementation of ctors and = operator are added
88 * Revision 1.96 2006/08/25 16:56:30 kharlov
89 * Compliance with Effective C++
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)
94 * Revision 1.94 2006/08/07 12:27:49 hristov
95 * Removing obsolete code which affected the event numbering scheme
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
100 * Revision 1.92 2006/04/29 20:26:46 hristov
101 * Separate EMC and CPV calibration (Yu.Kharlov)
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)
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)
109 * Revision 1.89 2006/03/13 14:05:42 kharlov
110 * Calibration objects for EMC and CPV
112 * Revision 1.88 2006/01/11 08:54:52 hristov
113 * Additional protection in case no calibration entry was found
115 * Revision 1.87 2005/11/22 08:46:43 kharlov
116 * Updated to new CDB (Boris Polichtchouk)
118 * Revision 1.86 2005/11/14 21:52:43 hristov
121 * Revision 1.85 2005/09/27 16:08:08 hristov
122 * New version of CDB storage framework (A.Colla)
124 * Revision 1.84 2005/09/21 10:02:47 kharlov
125 * Reading calibration from CDB (Boris Polichtchouk)
127 * Revision 1.82 2005/09/02 15:43:13 kharlov
128 * Add comments in GetCalibrationParameters and Calibrate
130 * Revision 1.81 2005/09/02 14:32:07 kharlov
131 * Calibration of raw data
133 * Revision 1.80 2005/08/24 15:31:36 kharlov
134 * Setting raw digits flag
136 * Revision 1.79 2005/07/25 15:53:53 kharlov
139 * Revision 1.78 2005/05/28 14:19:04 schutz
140 * Compilation warnings fixed by T.P.
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
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
164 // --- ROOT system ---
169 #include "TBenchmark.h"
170 #include "TClonesArray.h"
172 // --- Standard library ---
174 // --- AliRoot header files ---
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 "AliCDBManager.h"
184 #include "AliCDBStorage.h"
185 #include "AliCDBEntry.h"
186 #include "AliPHOSRecoParam.h"
187 #include "AliPHOSReconstructor.h"
188 #include "AliPHOSCalibData.h"
190 ClassImp(AliPHOSClusterizerv1)
192 //____________________________________________________________________________
193 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
194 AliPHOSClusterizer(),
195 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
197 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
198 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
199 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
200 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
202 // default ctor (to be used mainly by Streamer)
204 fDefaultInit = kTRUE ;
207 //____________________________________________________________________________
208 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
209 AliPHOSClusterizer(geom),
210 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
212 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
213 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
214 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
215 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
217 // ctor with the indication of the file where header Tree and digits Tree are stored
220 fDefaultInit = kFALSE ;
223 //____________________________________________________________________________
224 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
229 //____________________________________________________________________________
230 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
232 // Steering method to perform clusterization for one event
233 // The input is the tree with digits.
234 // The output is the tree with clusters.
236 if(strstr(option,"tim"))
237 gBenchmark->Start("PHOSClusterizer");
239 if(strstr(option,"print")) {
246 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
247 fEMCRecPoints->GetEntries()));
248 if(AliLog::GetGlobalDebugLevel()>1)
249 fEMCRecPoints->Print();
256 if(strstr(option,"deb"))
257 PrintRecPoints(option) ;
259 if(strstr(option,"tim")){
260 gBenchmark->Stop("PHOSClusterizer");
261 AliInfo(Form("took %f seconds for Clusterizing\n",
262 gBenchmark->GetCpuTime("PHOSClusterizer")));
264 fEMCRecPoints->Delete();
265 fCPVRecPoints->Delete();
268 //____________________________________________________________________________
269 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
270 Int_t nPar, Float_t * fitparameters) const
272 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
273 // The initial values for fitting procedure are set equal to the positions of local maxima.
274 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
277 gMinuit->mncler(); // Reset Minuit's list of paramters
278 gMinuit->SetPrintLevel(-1) ; // No Printout
279 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
280 // To set the address of the minimization function
282 TList * toMinuit = new TList();
283 toMinuit->AddAt(emcRP,0) ;
284 toMinuit->AddAt(fDigitsArr,1) ;
285 toMinuit->AddAt(fGeom,2) ;
287 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
289 // filling initial values for fit parameters
290 AliPHOSDigit * digit ;
294 Int_t nDigits = (Int_t) nPar / 3 ;
298 for(iDigit = 0; iDigit < nDigits; iDigit++){
299 digit = maxAt[iDigit];
304 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
305 fGeom->RelPosInModule(relid, x, z) ;
307 Float_t energy = maxAtEnergy[iDigit] ;
309 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
312 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
315 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
318 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
321 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
324 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
329 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
334 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
335 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
336 gMinuit->SetMaxIterations(5);
337 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
339 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
341 if(ierflg == 4){ // Minimum not found
342 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
345 for(index = 0; index < nPar; index++){
348 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
349 fitparameters[index] = val ;
358 //____________________________________________________________________________
359 void AliPHOSClusterizerv1::Init()
361 // Make all memory allocations which can not be done in default constructor.
362 // Attach the Clusterizer task to the list of PHOS tasks
364 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
367 gMinuit = new TMinuit(100);
370 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
371 if (fgCalibData->GetCalibDataEmc() == 0)
372 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
373 if (fgCalibData->GetCalibDataCpv() == 0)
374 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
378 //____________________________________________________________________________
379 void AliPHOSClusterizerv1::InitParameters()
382 fNumberOfCpvClusters = 0 ;
383 fNumberOfEmcClusters = 0 ;
385 const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
386 if(!recoParam) AliFatal("Reconstruction parameters are not set!");
390 fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
391 fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
393 fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
394 fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
396 fW0 = recoParam->GetEMCLogWeight();
397 fW0CPV = recoParam->GetCPVLogWeight();
399 fEmcTimeGate = 1.e-6 ; //10 sample steps
400 fEcoreRadius = recoParam->GetEMCEcoreRadius();
402 fToUnfold = recoParam->EMCToUnfold() ;
407 //____________________________________________________________________________
408 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
410 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
412 // = 2 are not neighbour but do not continue searching
413 // =-1 are not neighbour, continue searching, but do not look before d2 next time
414 // neighbours are defined as digits having at least a common vertex
415 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
416 // which is compared to a digit (d2) not yet in a cluster
419 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
422 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
424 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
425 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
426 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
428 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
429 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
430 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
434 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
435 return 2; // Difference in row numbers is too large to look further
441 if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
443 if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
452 //____________________________________________________________________________
453 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
455 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
459 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
461 if(digit->GetId() <= nEMC ) rv = kTRUE;
466 //____________________________________________________________________________
467 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
469 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
473 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
475 if(digit->GetId() > nEMC ) rv = kTRUE;
480 //____________________________________________________________________________
481 void AliPHOSClusterizerv1::WriteRecPoints()
484 // Creates new branches with given title
485 // fills and writes into TreeR.
488 //Evaluate position, dispersion and other RecPoint properties..
489 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
490 Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
491 TVector3 fakeVtx(0.,0.,0.) ;
492 for(index = 0; index < nEmc; index++){
493 AliPHOSEmcRecPoint * rp =
494 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
495 rp->Purify(emcMinE) ;
496 if(rp->GetMultiplicity()==0){
497 fEMCRecPoints->RemoveAt(index) ;
502 // No vertex is available now, calculate corrections in PID
503 rp->EvalAll(fDigitsArr) ;
504 rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
505 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
506 rp->EvalLocal2TrackingCSTransform();
508 fEMCRecPoints->Compress() ;
509 fEMCRecPoints->Sort() ;
510 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
511 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
512 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
515 //For each rec.point set the distance to the nearest bad crystal (BVP)
516 SetDistancesToBadChannels();
518 //Now the same for CPV
519 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
520 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
521 rp->EvalAll(fDigitsArr) ;
522 rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
523 rp->EvalLocal2TrackingCSTransform();
525 fCPVRecPoints->Sort() ;
527 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
528 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
530 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
532 if(fWrite){ //We write TreeR
537 //____________________________________________________________________________
538 void AliPHOSClusterizerv1::MakeClusters()
540 // Steering method to construct the clusters stored in a list of Reconstructed Points
541 // A cluster is defined as a list of neighbour digits
543 fNumberOfCpvClusters = 0 ;
544 fNumberOfEmcClusters = 0 ;
546 //Mark all digits as unused yet
547 const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
548 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
550 for(Int_t i=0; i<nDigits; i++){
553 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
554 //e.g. first digit in this module, first CPV digit etc.
555 AliPHOSDigit * digit ;
556 TArrayI clusterdigitslist(maxNDigits) ;
557 AliPHOSRecPoint * clu = 0 ;
558 for(Int_t i=0; i<nDigits; i++){
562 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
568 //is this digit so energetic that start cluster?
569 if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
570 ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
571 Int_t iDigitInCluster = 0 ;
572 if ( IsInEmc(digit) ) {
573 // start a new EMC RecPoint
574 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
575 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
577 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
578 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
579 fNumberOfEmcClusters++ ;
580 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId())) ;
581 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
583 fDigitsUsed[i]=kTRUE ;
585 // start a new CPV cluster
586 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
587 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
589 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
590 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
591 fNumberOfCpvClusters++ ;
592 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
593 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
595 fDigitsUsed[i]=kTRUE ;
598 //Now scan remaining digits in list to find neigbours of our seed
600 AliPHOSDigit * digitN ;
602 while (index < iDigitInCluster){ // scan over digits already in cluster
603 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
605 for(Int_t j=iFirst; j<nDigits; j++){
606 if (iDigitInCluster >= maxNDigits) {
607 AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
612 continue ; //look through remaining digits
613 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
614 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
616 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
619 case 0 : // not a neighbour
621 case 1 : // are neighbours
622 clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digit->GetId()),CalibrateT(digitN->GetTime(),digit->GetId())) ;
623 clusterdigitslist[iDigitInCluster] = j ;
625 fDigitsUsed[j]=kTRUE ;
627 case 2 : // too far from each other
633 endOfLoop: ; //scanned all possible neighbours for this digit
635 } // loop over cluster
641 //____________________________________________________________________________
642 void AliPHOSClusterizerv1::MakeUnfolding()
644 // Unfolds clusters using the shape of an ElectroMagnetic shower
645 // Performs unfolding of all EMC/CPV clusters
647 // Unfold first EMC clusters
648 if(fNumberOfEmcClusters > 0){
650 Int_t nModulesToUnfold = fGeom->GetNModules() ;
652 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
654 for(index = 0 ; index < numberofNotUnfolded ; index++){
656 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
657 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
660 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
661 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
662 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
663 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
665 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
666 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
668 fEMCRecPoints->Remove(emcRecPoint);
669 fEMCRecPoints->Compress() ;
671 fNumberOfEmcClusters -- ;
672 numberofNotUnfolded-- ;
675 emcRecPoint->SetNExMax(1) ; //Only one local maximum
679 delete[] maxAtEnergy ;
682 // Unfolding of EMC clusters finished
685 // Unfold now CPV clusters
686 if(fNumberOfCpvClusters > 0){
688 Int_t nModulesToUnfold = fGeom->GetNModules() ;
690 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
692 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
694 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
696 if(recPoint->GetPHOSMod()> nModulesToUnfold)
699 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
701 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
702 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
703 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
704 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
706 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
707 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
708 fCPVRecPoints->Remove(emcRecPoint);
709 fCPVRecPoints->Compress() ;
711 numberofCpvNotUnfolded-- ;
712 fNumberOfCpvClusters-- ;
716 delete[] maxAtEnergy ;
719 //Unfolding of Cpv clusters finished
723 //____________________________________________________________________________
724 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
726 // Shape of the shower (see PHOS TDR)
727 // If you change this function, change also the gradient evaluation in ChiSquare()
729 //for the moment we neglect dependence on the incident angle.
731 Double_t r2 = x*x + z*z ;
732 Double_t r4 = r2*r2 ;
733 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
734 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
738 //____________________________________________________________________________
739 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
741 AliPHOSDigit ** maxAt,
742 Float_t * maxAtEnergy)
744 // Performs the unfolding of a cluster with nMax overlapping showers
746 Int_t nPar = 3 * nMax ;
747 Float_t * fitparameters = new Float_t[nPar] ;
749 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
752 // Fit failed, return and remove cluster
753 iniEmc->SetNExMax(-1) ;
754 delete[] fitparameters ;
758 // create ufolded rec points and fill them with new energy lists
759 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
760 // and later correct this number in acordance with actual energy deposition
762 Int_t nDigits = iniEmc->GetMultiplicity() ;
763 Float_t * efit = new Float_t[nDigits] ;
764 Float_t xDigit=0.,zDigit=0. ;
765 Float_t xpar=0.,zpar=0.,epar=0. ;
767 AliPHOSDigit * digit = 0 ;
768 Int_t * emcDigits = iniEmc->GetDigitsList() ;
774 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
775 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
776 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
777 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
781 while(iparam < nPar ){
782 xpar = fitparameters[iparam] ;
783 zpar = fitparameters[iparam+1] ;
784 epar = fitparameters[iparam+2] ;
786 // fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
787 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
788 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
792 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
793 // so that energy deposited in each cell is distributed betwin new clusters proportionally
794 // to its contribution to efit
796 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
800 while(iparam < nPar ){
801 xpar = fitparameters[iparam] ;
802 zpar = fitparameters[iparam+1] ;
803 epar = fitparameters[iparam+2] ;
805 // fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
807 AliPHOSEmcRecPoint * emcRP = 0 ;
809 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
811 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
812 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
814 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
815 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
816 fNumberOfEmcClusters++ ;
817 emcRP->SetNExMax((Int_t)nPar/3) ;
819 else{//create new entries in fCPVRecPoints
820 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
821 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
823 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
824 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
825 fNumberOfCpvClusters++ ;
829 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
830 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
831 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
832 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
833 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
834 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
835 eDigit = emcEnergies[iDigit] * ratio ;
836 emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId()) ) ;
840 delete[] fitparameters ;
845 //_____________________________________________________________________________
846 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
848 // Calculates the Chi square for the cluster unfolding minimization
849 // Number of parameters, Gradient, Chi squared, parameters, what to do
851 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
853 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
854 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
855 // A bit buggy way to get an access to the geometry
857 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
859 // TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
861 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
863 Int_t * emcDigits = emcRP->GetDigitsList() ;
865 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
867 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
874 for(iparam = 0 ; iparam < nPar ; iparam++)
875 Grad[iparam] = 0 ; // Will evaluate gradient
879 AliPHOSDigit * digit ;
882 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
884 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
890 geom->AbsToRelNumbering(digit->GetId(), relid) ;
892 geom->RelPosInModule(relid, xDigit, zDigit) ;
894 if(iflag == 2){ // calculate gradient
897 while(iParam < nPar ){
898 Double_t dx = (xDigit - x[iParam]) ;
900 Double_t dz = (zDigit - x[iParam]) ;
902 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
903 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
904 efit += x[iParam] * ShowerShape(dx,dz) ;
907 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
909 while(iParam < nPar ){
910 Double_t xpar = x[iParam] ;
911 Double_t zpar = x[iParam+1] ;
912 Double_t epar = x[iParam+2] ;
913 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
914 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
915 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
916 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
917 //DP: No incident angle dependence in gradient yet!!!!!!
918 Double_t r4 = dr*dr*dr*dr ;
919 Double_t r295 = TMath::Power(dr,2.95) ;
920 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
921 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
923 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
925 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
927 Grad[iParam] += shape ; // Derivative over energy
934 while(iparam < nPar ){
935 Double_t xpar = x[iparam] ;
936 Double_t zpar = x[iparam+1] ;
937 Double_t epar = x[iparam+2] ;
939 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
940 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
941 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
944 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
945 // Here we assume, that sigma = sqrt(E)
950 //____________________________________________________________________________
951 void AliPHOSClusterizerv1::Print(const Option_t *)const
953 // Print clusterizer parameters
956 TString taskName(GetName()) ;
957 taskName.ReplaceAll(Version(), "") ;
959 if( strcmp(GetName(), "") !=0 ) {
961 message = "\n--------------- %s %s -----------\n" ;
962 message += "Clusterizing digits from the file: %s\n" ;
963 message += " Branch: %s\n" ;
964 message += " EMC Clustering threshold = %f\n" ;
965 message += " EMC Local Maximum cut = %f\n" ;
966 message += " EMC Logarothmic weight = %f\n" ;
967 message += " CPV Clustering threshold = %f\n" ;
968 message += " CPV Local Maximum cut = %f\n" ;
969 message += " CPV Logarothmic weight = %f\n" ;
971 message += " Unfolding on\n" ;
973 message += " Unfolding off\n" ;
975 message += "------------------------------------------------------------------" ;
978 message = " AliPHOSClusterizerv1 not initialized " ;
980 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
985 fEmcClusteringThreshold,
988 fCpvClusteringThreshold,
992 //____________________________________________________________________________
993 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
995 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
997 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
998 fEMCRecPoints->GetEntriesFast(),
999 fCPVRecPoints->GetEntriesFast() )) ;
1001 if(strstr(option,"all")) {
1002 printf("\n EMC clusters \n") ;
1003 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1005 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1006 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1008 rp->GetLocalPosition(locpos);
1010 rp->GetElipsAxis(lambda);
1013 primaries = rp->GetPrimaries(nprimaries);
1014 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1015 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1016 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1018 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1019 printf("%d ", primaries[iprimary] ) ;
1024 //Now plot CPV recPoints
1025 printf("\n CPV clusters \n") ;
1026 printf("Index Ene(MeV) Module X Y Z \n") ;
1027 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1028 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1031 rp->GetLocalPosition(locpos);
1033 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1034 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1035 locpos.X(), locpos.Y(), locpos.Z()) ;
1041 //____________________________________________________________________________
1042 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1044 //For each EMC rec. point set the distance to the nearest bad crystal.
1045 //Author: Boris Polichtchouk
1047 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1048 AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1051 fgCalibData->EmcBadChannelIds(badIds);
1053 AliPHOSEmcRecPoint* rp;
1056 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1057 TVector3 gposBadChannel; // global position of bad crystal
1060 Float_t dist,minDist;
1061 Int_t relid[4]={0,0,0,0} ;
1063 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1064 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1065 //evaluate distance to border
1066 relid[0]=rp->GetPHOSMod() ;
1069 Float_t xcorner,zcorner;
1070 fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1071 rp->GetLocalPosition(lpos) ;
1072 minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1073 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1074 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1075 if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1076 continue ; //bad channels can be in the module which does not exist in simulations.
1077 rp->GetGlobalPosition(gposRecPoint,gmat);
1078 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1079 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1080 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1081 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1082 dR = gposBadChannel-gposRecPoint;
1084 if(dist<minDist) minDist = dist;
1087 rp->SetDistanceToBadCrystal(minDist);
1091 //==================================================================================
1092 Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
1093 // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1095 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1097 //Determine rel.position of the cell absolute ID
1099 geom->AbsToRelNumbering(absId,relId);
1100 Int_t module=relId[0];
1101 Int_t row =relId[2];
1102 Int_t column=relId[3];
1104 Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1105 return amp*calibration ;
1108 Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1109 return amp*calibration ;
1112 //==================================================================================
1113 Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId)const{
1114 // Calibrate time in EMC digit
1116 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1118 //Determine rel.position of the cell absolute ID
1120 geom->AbsToRelNumbering(absId,relId);
1121 Int_t module=relId[0];
1122 Int_t row =relId[2];
1123 Int_t column=relId[3];
1128 time += fgCalibData->GetTimeShiftEmc(module,column,row);