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 "AliPHOSCalibrationDB.h"
184 #include "AliCDBManager.h"
185 #include "AliCDBStorage.h"
186 #include "AliCDBEntry.h"
187 #include "AliPHOSRecoParam.h"
188 #include "AliPHOSReconstructor.h"
189 #include "AliPHOSCalibData.h"
191 ClassImp(AliPHOSClusterizerv1)
193 //____________________________________________________________________________
194 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
195 AliPHOSClusterizer(),
196 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
198 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
199 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
200 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
201 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
203 // default ctor (to be used mainly by Streamer)
205 fDefaultInit = kTRUE ;
208 //____________________________________________________________________________
209 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
210 AliPHOSClusterizer(geom),
211 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
213 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
214 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
215 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
216 fW0CPV(0), fEmcTimeGate(0), fEcoreRadius(0)
218 // ctor with the indication of the file where header Tree and digits Tree are stored
221 fDefaultInit = kFALSE ;
224 //____________________________________________________________________________
225 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
230 //____________________________________________________________________________
231 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
233 // Steering method to perform clusterization for one event
234 // The input is the tree with digits.
235 // The output is the tree with clusters.
237 if(strstr(option,"tim"))
238 gBenchmark->Start("PHOSClusterizer");
240 if(strstr(option,"print")) {
247 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
248 fEMCRecPoints->GetEntries()));
249 if(AliLog::GetGlobalDebugLevel()>1)
250 fEMCRecPoints->Print();
257 if(strstr(option,"deb"))
258 PrintRecPoints(option) ;
260 if(strstr(option,"tim")){
261 gBenchmark->Stop("PHOSClusterizer");
262 AliInfo(Form("took %f seconds for Clusterizing\n",
263 gBenchmark->GetCpuTime("PHOSClusterizer")));
265 fEMCRecPoints->Delete();
266 fCPVRecPoints->Delete();
269 //____________________________________________________________________________
270 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
271 Int_t nPar, Float_t * fitparameters) const
273 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
274 // The initial values for fitting procedure are set equal to the positions of local maxima.
275 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
278 gMinuit->mncler(); // Reset Minuit's list of paramters
279 gMinuit->SetPrintLevel(-1) ; // No Printout
280 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
281 // To set the address of the minimization function
283 TList * toMinuit = new TList();
284 toMinuit->AddAt(emcRP,0) ;
285 toMinuit->AddAt(fDigitsArr,1) ;
286 toMinuit->AddAt(fGeom,2) ;
288 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
290 // filling initial values for fit parameters
291 AliPHOSDigit * digit ;
295 Int_t nDigits = (Int_t) nPar / 3 ;
299 for(iDigit = 0; iDigit < nDigits; iDigit++){
300 digit = maxAt[iDigit];
305 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
306 fGeom->RelPosInModule(relid, x, z) ;
308 Float_t energy = maxAtEnergy[iDigit] ;
310 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
313 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
316 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
319 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
322 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
325 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
330 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
335 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
336 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
337 gMinuit->SetMaxIterations(5);
338 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
340 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
342 if(ierflg == 4){ // Minimum not found
343 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
346 for(index = 0; index < nPar; index++){
349 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
350 fitparameters[index] = val ;
359 //____________________________________________________________________________
360 void AliPHOSClusterizerv1::Init()
362 // Make all memory allocations which can not be done in default constructor.
363 // Attach the Clusterizer task to the list of PHOS tasks
365 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
368 gMinuit = new TMinuit(100);
371 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
372 if (fgCalibData->GetCalibDataEmc() == 0)
373 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
377 //____________________________________________________________________________
378 void AliPHOSClusterizerv1::InitParameters()
381 fNumberOfCpvClusters = 0 ;
382 fNumberOfEmcClusters = 0 ;
384 const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
385 if(!recoParam) AliFatal("Reconstruction parameters are not set!");
389 fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
390 fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
392 fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
393 fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
395 fW0 = recoParam->GetEMCLogWeight();
396 fW0CPV = recoParam->GetCPVLogWeight();
398 fEmcTimeGate = 1.e-6 ;
399 fEcoreRadius = recoParam->GetEMCEcoreRadius();
401 fToUnfold = recoParam->EMCToUnfold() ;
406 //____________________________________________________________________________
407 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
409 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
411 // = 2 are not neighbour but do not continue searching
412 // =-1 are not neighbour, continue searching, but do not look before d2 next time
413 // neighbours are defined as digits having at least a common vertex
414 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
415 // which is compared to a digit (d2) not yet in a cluster
418 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
421 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
423 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
424 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
425 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
427 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
428 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
429 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
433 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
434 return 2; // Difference in row numbers is too large to look further
440 if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
442 if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
451 //____________________________________________________________________________
452 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
454 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
458 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
460 if(digit->GetId() <= nEMC ) rv = kTRUE;
465 //____________________________________________________________________________
466 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
468 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
472 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
474 if(digit->GetId() > nEMC ) rv = kTRUE;
479 //____________________________________________________________________________
480 void AliPHOSClusterizerv1::WriteRecPoints()
483 // Creates new branches with given title
484 // fills and writes into TreeR.
487 //Evaluate position, dispersion and other RecPoint properties..
488 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
489 Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
490 TVector3 fakeVtx(0.,0.,0.) ;
491 for(index = 0; index < nEmc; index++){
492 AliPHOSEmcRecPoint * rp =
493 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
494 rp->Purify(emcMinE) ;
495 if(rp->GetMultiplicity()==0){
496 fEMCRecPoints->RemoveAt(index) ;
501 // No vertex is available now, calculate corrections in PID
502 rp->EvalAll(fDigitsArr) ;
503 rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
504 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
505 rp->EvalLocal2TrackingCSTransform();
507 fEMCRecPoints->Compress() ;
508 fEMCRecPoints->Sort() ;
509 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
510 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
511 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
514 //For each rec.point set the distance to the nearest bad crystal (BVP)
515 SetDistancesToBadChannels();
517 //Now the same for CPV
518 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
519 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
520 rp->EvalAll(fDigitsArr) ;
521 rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
522 rp->EvalLocal2TrackingCSTransform();
524 fCPVRecPoints->Sort() ;
526 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
527 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
529 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
531 if(fWrite){ //We write TreeR
536 //____________________________________________________________________________
537 void AliPHOSClusterizerv1::MakeClusters()
539 // Steering method to construct the clusters stored in a list of Reconstructed Points
540 // A cluster is defined as a list of neighbour digits
542 fNumberOfCpvClusters = 0 ;
543 fNumberOfEmcClusters = 0 ;
545 //Mark all digits as unused yet
546 const Int_t maxNDigits = 1500;
547 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
548 if (nDigits > maxNDigits) {
549 AliWarning(Form("Skip event with too high digit occupancy: nDigits=%d",nDigits));
553 for(Int_t i=0; i<nDigits; i++){
556 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
557 //e.g. first digit in this module, first CPV digit etc.
558 AliPHOSDigit * digit ;
559 TArrayI clusterdigitslist(maxNDigits) ;
560 AliPHOSRecPoint * clu = 0 ;
561 for(Int_t i=0; i<nDigits; i++){
565 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
571 //is this digit so energetic that start cluster?
572 if (( IsInEmc(digit) && digit->GetEnergy() > fEmcClusteringThreshold ) ||
573 ( IsInCpv(digit) && digit->GetEnergy() > fCpvClusteringThreshold ) ) {
574 Int_t iDigitInCluster = 0 ;
576 if ( IsInEmc(digit) ) {
577 // start a new EMC RecPoint
578 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
579 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
581 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
582 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
583 fNumberOfEmcClusters++ ;
584 clu->AddDigit(*digit, digit->GetEnergy()) ;
585 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
587 fDigitsUsed[i]=kTRUE ;
589 // start a new CPV cluster
590 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
591 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
593 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
594 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
595 fNumberOfCpvClusters++ ;
596 clu->AddDigit(*digit, digit->GetEnergy()) ;
597 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
599 fDigitsUsed[i]=kTRUE ;
602 //Now scan remaining digits in list to find neigbours of our seed
604 AliPHOSDigit * digitN ;
606 while (index < iDigitInCluster){ // scan over digits already in cluster
607 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
609 for(Int_t j=iFirst; j<nDigits; j++){
611 continue ; //look through remaining digits
612 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
613 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
615 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
618 case 0 : // not a neighbour
620 case 1 : // are neighbours
621 clu->AddDigit(*digitN, digitN->GetEnergy());
622 clusterdigitslist[iDigitInCluster] = j ;
624 fDigitsUsed[j]=kTRUE ;
626 case 2 : // too far from each other
632 endofloop: ; //scanned all possible neighbours for this digit
634 } // loop over cluster
640 //____________________________________________________________________________
641 void AliPHOSClusterizerv1::MakeUnfolding()
643 // Unfolds clusters using the shape of an ElectroMagnetic shower
644 // Performs unfolding of all EMC/CPV clusters
646 // Unfold first EMC clusters
647 if(fNumberOfEmcClusters > 0){
649 Int_t nModulesToUnfold = fGeom->GetNModules() ;
651 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
653 for(index = 0 ; index < numberofNotUnfolded ; index++){
655 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
656 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
659 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
660 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
661 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
662 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
664 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
665 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
667 fEMCRecPoints->Remove(emcRecPoint);
668 fEMCRecPoints->Compress() ;
670 fNumberOfEmcClusters -- ;
671 numberofNotUnfolded-- ;
674 emcRecPoint->SetNExMax(1) ; //Only one local maximum
678 delete[] maxAtEnergy ;
681 // Unfolding of EMC clusters finished
684 // Unfold now CPV clusters
685 if(fNumberOfCpvClusters > 0){
687 Int_t nModulesToUnfold = fGeom->GetNModules() ;
689 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
691 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
693 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
695 if(recPoint->GetPHOSMod()> nModulesToUnfold)
698 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
700 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
701 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
702 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
703 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
705 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
706 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
707 fCPVRecPoints->Remove(emcRecPoint);
708 fCPVRecPoints->Compress() ;
710 numberofCpvNotUnfolded-- ;
711 fNumberOfCpvClusters-- ;
715 delete[] maxAtEnergy ;
718 //Unfolding of Cpv clusters finished
722 //____________________________________________________________________________
723 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
725 // Shape of the shower (see PHOS TDR)
726 // If you change this function, change also the gradient evaluation in ChiSquare()
728 //for the moment we neglect dependence on the incident angle.
730 Double_t r2 = x*x + z*z ;
731 Double_t r4 = r2*r2 ;
732 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
733 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
737 //____________________________________________________________________________
738 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
740 AliPHOSDigit ** maxAt,
741 Float_t * maxAtEnergy)
743 // Performs the unfolding of a cluster with nMax overlapping showers
745 Int_t nPar = 3 * nMax ;
746 Float_t * fitparameters = new Float_t[nPar] ;
748 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
751 // Fit failed, return and remove cluster
752 iniEmc->SetNExMax(-1) ;
753 delete[] fitparameters ;
757 // create ufolded rec points and fill them with new energy lists
758 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
759 // and later correct this number in acordance with actual energy deposition
761 Int_t nDigits = iniEmc->GetMultiplicity() ;
762 Float_t * efit = new Float_t[nDigits] ;
763 Float_t xDigit=0.,zDigit=0. ;
764 Float_t xpar=0.,zpar=0.,epar=0. ;
766 AliPHOSDigit * digit = 0 ;
767 Int_t * emcDigits = iniEmc->GetDigitsList() ;
773 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
774 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
775 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
776 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
780 while(iparam < nPar ){
781 xpar = fitparameters[iparam] ;
782 zpar = fitparameters[iparam+1] ;
783 epar = fitparameters[iparam+2] ;
785 // fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
786 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
787 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
791 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
792 // so that energy deposited in each cell is distributed betwin new clusters proportionally
793 // to its contribution to efit
795 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
799 while(iparam < nPar ){
800 xpar = fitparameters[iparam] ;
801 zpar = fitparameters[iparam+1] ;
802 epar = fitparameters[iparam+2] ;
804 // fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
806 AliPHOSEmcRecPoint * emcRP = 0 ;
808 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
810 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
811 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
813 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
814 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
815 fNumberOfEmcClusters++ ;
816 emcRP->SetNExMax((Int_t)nPar/3) ;
818 else{//create new entries in fCPVRecPoints
819 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
820 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
822 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
823 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
824 fNumberOfCpvClusters++ ;
828 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
829 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
830 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
831 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
832 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
833 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
834 eDigit = emcEnergies[iDigit] * ratio ;
835 emcRP->AddDigit( *digit, eDigit ) ;
839 delete[] fitparameters ;
844 //_____________________________________________________________________________
845 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
847 // Calculates the Chi square for the cluster unfolding minimization
848 // Number of parameters, Gradient, Chi squared, parameters, what to do
850 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
852 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
853 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
854 // A bit buggy way to get an access to the geometry
856 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
858 // TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
860 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
862 Int_t * emcDigits = emcRP->GetDigitsList() ;
864 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
866 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
873 for(iparam = 0 ; iparam < nPar ; iparam++)
874 Grad[iparam] = 0 ; // Will evaluate gradient
878 AliPHOSDigit * digit ;
881 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
883 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
889 geom->AbsToRelNumbering(digit->GetId(), relid) ;
891 geom->RelPosInModule(relid, xDigit, zDigit) ;
893 if(iflag == 2){ // calculate gradient
896 while(iParam < nPar ){
897 Double_t dx = (xDigit - x[iParam]) ;
899 Double_t dz = (zDigit - x[iParam]) ;
901 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
902 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
903 efit += x[iParam] * ShowerShape(dx,dz) ;
906 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
908 while(iParam < nPar ){
909 Double_t xpar = x[iParam] ;
910 Double_t zpar = x[iParam+1] ;
911 Double_t epar = x[iParam+2] ;
912 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
913 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
914 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
915 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
916 //DP: No incident angle dependence in gradient yet!!!!!!
917 Double_t r4 = dr*dr*dr*dr ;
918 Double_t r295 = TMath::Power(dr,2.95) ;
919 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
920 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
922 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
924 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
926 Grad[iParam] += shape ; // Derivative over energy
933 while(iparam < nPar ){
934 Double_t xpar = x[iparam] ;
935 Double_t zpar = x[iparam+1] ;
936 Double_t epar = x[iparam+2] ;
938 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
939 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
940 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
943 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
944 // Here we assume, that sigma = sqrt(E)
949 //____________________________________________________________________________
950 void AliPHOSClusterizerv1::Print(const Option_t *)const
952 // Print clusterizer parameters
955 TString taskName(GetName()) ;
956 taskName.ReplaceAll(Version(), "") ;
958 if( strcmp(GetName(), "") !=0 ) {
960 message = "\n--------------- %s %s -----------\n" ;
961 message += "Clusterizing digits from the file: %s\n" ;
962 message += " Branch: %s\n" ;
963 message += " EMC Clustering threshold = %f\n" ;
964 message += " EMC Local Maximum cut = %f\n" ;
965 message += " EMC Logarothmic weight = %f\n" ;
966 message += " CPV Clustering threshold = %f\n" ;
967 message += " CPV Local Maximum cut = %f\n" ;
968 message += " CPV Logarothmic weight = %f\n" ;
970 message += " Unfolding on\n" ;
972 message += " Unfolding off\n" ;
974 message += "------------------------------------------------------------------" ;
977 message = " AliPHOSClusterizerv1 not initialized " ;
979 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
984 fEmcClusteringThreshold,
987 fCpvClusteringThreshold,
991 //____________________________________________________________________________
992 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
994 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
996 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
997 fEMCRecPoints->GetEntriesFast(),
998 fCPVRecPoints->GetEntriesFast() )) ;
1000 if(strstr(option,"all")) {
1001 printf("\n EMC clusters \n") ;
1002 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1004 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1005 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1007 rp->GetLocalPosition(locpos);
1009 rp->GetElipsAxis(lambda);
1012 primaries = rp->GetPrimaries(nprimaries);
1013 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1014 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1015 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1017 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1018 printf("%d ", primaries[iprimary] ) ;
1023 //Now plot CPV recPoints
1024 printf("\n CPV clusters \n") ;
1025 printf("Index Ene(MeV) Module X Y Z \n") ;
1026 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1027 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1030 rp->GetLocalPosition(locpos);
1032 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1033 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1034 locpos.X(), locpos.Y(), locpos.Z()) ;
1040 //____________________________________________________________________________
1041 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1043 //For each EMC rec. point set the distance to the nearest bad crystal.
1044 //Author: Boris Polichtchouk
1046 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1047 AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1050 fgCalibData->EmcBadChannelIds(badIds);
1052 AliPHOSEmcRecPoint* rp;
1055 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1056 TVector3 gposBadChannel; // global position of bad crystal
1059 Float_t dist,minDist;
1060 Int_t relid[4]={0,0,0,0} ;
1062 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1063 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1064 //evaluate distance to border
1065 relid[0]=rp->GetPHOSMod() ;
1068 Float_t xcorner,zcorner;
1069 fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1070 rp->GetLocalPosition(lpos) ;
1071 minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1072 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1073 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1074 if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1075 continue ; //bad channels can be in the module which does not exist in simulations.
1076 rp->GetGlobalPosition(gposRecPoint,gmat);
1077 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1078 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1079 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1080 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1081 dR = gposBadChannel-gposRecPoint;
1083 if(dist<minDist) minDist = dist;
1086 rp->SetDistanceToBadCrystal(minDist);