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 if(!gMinuit) //it was deleted by someone else
278 gMinuit = new TMinuit(100) ;
279 gMinuit->mncler(); // Reset Minuit's list of paramters
280 gMinuit->SetPrintLevel(-1) ; // No Printout
281 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
282 // To set the address of the minimization function
284 TList * toMinuit = new TList();
285 toMinuit->AddAt(emcRP,0) ;
286 toMinuit->AddAt(fDigitsArr,1) ;
287 toMinuit->AddAt(fGeom,2) ;
289 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
291 // filling initial values for fit parameters
292 AliPHOSDigit * digit ;
296 Int_t nDigits = (Int_t) nPar / 3 ;
300 for(iDigit = 0; iDigit < nDigits; iDigit++){
301 digit = maxAt[iDigit];
306 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
307 fGeom->RelPosInModule(relid, x, z) ;
309 Float_t energy = maxAtEnergy[iDigit] ;
311 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
314 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
317 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
320 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
323 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
326 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
331 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
336 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
337 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
338 gMinuit->SetMaxIterations(5);
339 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
341 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
343 if(ierflg == 4){ // Minimum not found
344 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
347 for(index = 0; index < nPar; index++){
350 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
351 fitparameters[index] = val ;
360 //____________________________________________________________________________
361 void AliPHOSClusterizerv1::Init()
363 // Make all memory allocations which can not be done in default constructor.
364 // Attach the Clusterizer task to the list of PHOS tasks
366 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
369 gMinuit = new TMinuit(100);
372 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
373 if (fgCalibData->GetCalibDataEmc() == 0)
374 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
375 if (fgCalibData->GetCalibDataCpv() == 0)
376 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
380 //____________________________________________________________________________
381 void AliPHOSClusterizerv1::InitParameters()
384 fNumberOfCpvClusters = 0 ;
385 fNumberOfEmcClusters = 0 ;
387 const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
388 if(!recoParam) AliFatal("Reconstruction parameters are not set!");
392 fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
393 fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
395 fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
396 fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
398 fW0 = recoParam->GetEMCLogWeight();
399 fW0CPV = recoParam->GetCPVLogWeight();
401 fEmcTimeGate = 1.e-6 ; //10 sample steps
402 fEcoreRadius = recoParam->GetEMCEcoreRadius();
404 fToUnfold = recoParam->EMCToUnfold() ;
409 //____________________________________________________________________________
410 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
412 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
414 // = 2 are not neighbour but do not continue searching
415 // =-1 are not neighbour, continue searching, but do not look before d2 next time
416 // neighbours are defined as digits having at least a common vertex
417 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
418 // which is compared to a digit (d2) not yet in a cluster
421 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
424 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
426 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
427 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
428 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
430 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
431 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
432 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
436 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
437 return 2; // Difference in row numbers is too large to look further
443 if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
445 if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
454 //____________________________________________________________________________
455 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
457 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
461 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
463 if(digit->GetId() <= nEMC ) rv = kTRUE;
468 //____________________________________________________________________________
469 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
471 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
475 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
477 if(digit->GetId() > nEMC ) rv = kTRUE;
482 //____________________________________________________________________________
483 void AliPHOSClusterizerv1::WriteRecPoints()
486 // Creates new branches with given title
487 // fills and writes into TreeR.
490 //Evaluate position, dispersion and other RecPoint properties..
491 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
492 Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
493 TVector3 fakeVtx(0.,0.,0.) ;
494 for(index = 0; index < nEmc; index++){
495 AliPHOSEmcRecPoint * rp =
496 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
497 rp->Purify(emcMinE) ;
498 if(rp->GetMultiplicity()==0){
499 fEMCRecPoints->RemoveAt(index) ;
504 // No vertex is available now, calculate corrections in PID
505 rp->EvalAll(fDigitsArr) ;
506 rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
507 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
508 rp->EvalLocal2TrackingCSTransform();
510 fEMCRecPoints->Compress() ;
511 fEMCRecPoints->Sort() ;
512 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
513 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
514 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
517 //For each rec.point set the distance to the nearest bad crystal (BVP)
518 SetDistancesToBadChannels();
520 //Now the same for CPV
521 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
522 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
523 rp->EvalAll(fDigitsArr) ;
524 rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
525 rp->EvalLocal2TrackingCSTransform();
527 fCPVRecPoints->Sort() ;
529 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
530 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
532 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
534 if(fWrite){ //We write TreeR
539 //____________________________________________________________________________
540 void AliPHOSClusterizerv1::MakeClusters()
542 // Steering method to construct the clusters stored in a list of Reconstructed Points
543 // A cluster is defined as a list of neighbour digits
545 fNumberOfCpvClusters = 0 ;
546 fNumberOfEmcClusters = 0 ;
548 //Mark all digits as unused yet
549 const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
550 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
552 for(Int_t i=0; i<nDigits; i++){
555 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
556 //e.g. first digit in this module, first CPV digit etc.
557 AliPHOSDigit * digit ;
558 TArrayI clusterdigitslist(maxNDigits) ;
559 AliPHOSRecPoint * clu = 0 ;
560 for(Int_t i=0; i<nDigits; i++){
564 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
570 //is this digit so energetic that start cluster?
571 if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
572 ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
573 Int_t iDigitInCluster = 0 ;
574 if ( IsInEmc(digit) ) {
575 // start a new EMC RecPoint
576 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
577 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
579 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
580 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
581 fNumberOfEmcClusters++ ;
582 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId())) ;
583 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
585 fDigitsUsed[i]=kTRUE ;
587 // start a new CPV cluster
588 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
589 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
591 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
592 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
593 fNumberOfCpvClusters++ ;
594 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
595 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
597 fDigitsUsed[i]=kTRUE ;
600 //Now scan remaining digits in list to find neigbours of our seed
602 AliPHOSDigit * digitN ;
604 while (index < iDigitInCluster){ // scan over digits already in cluster
605 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
607 for(Int_t j=iFirst; j<nDigits; j++){
608 if (iDigitInCluster >= maxNDigits) {
609 AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
614 continue ; //look through remaining digits
615 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
616 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
618 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
621 case 0 : // not a neighbour
623 case 1 : // are neighbours
624 clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digit->GetId()),CalibrateT(digitN->GetTime(),digit->GetId())) ;
625 clusterdigitslist[iDigitInCluster] = j ;
627 fDigitsUsed[j]=kTRUE ;
629 case 2 : // too far from each other
635 endOfLoop: ; //scanned all possible neighbours for this digit
637 } // loop over cluster
643 //____________________________________________________________________________
644 void AliPHOSClusterizerv1::MakeUnfolding()
646 // Unfolds clusters using the shape of an ElectroMagnetic shower
647 // Performs unfolding of all EMC/CPV clusters
649 // Unfold first EMC clusters
650 if(fNumberOfEmcClusters > 0){
652 Int_t nModulesToUnfold = fGeom->GetNModules() ;
654 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
656 for(index = 0 ; index < numberofNotUnfolded ; index++){
658 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
659 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
662 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
663 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
664 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
665 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
667 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
668 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
670 fEMCRecPoints->Remove(emcRecPoint);
671 fEMCRecPoints->Compress() ;
673 fNumberOfEmcClusters -- ;
674 numberofNotUnfolded-- ;
677 emcRecPoint->SetNExMax(1) ; //Only one local maximum
681 delete[] maxAtEnergy ;
684 // Unfolding of EMC clusters finished
687 // Unfold now CPV clusters
688 if(fNumberOfCpvClusters > 0){
690 Int_t nModulesToUnfold = fGeom->GetNModules() ;
692 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
694 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
696 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
698 if(recPoint->GetPHOSMod()> nModulesToUnfold)
701 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
703 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
704 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
705 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
706 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
708 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
709 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
710 fCPVRecPoints->Remove(emcRecPoint);
711 fCPVRecPoints->Compress() ;
713 numberofCpvNotUnfolded-- ;
714 fNumberOfCpvClusters-- ;
718 delete[] maxAtEnergy ;
721 //Unfolding of Cpv clusters finished
725 //____________________________________________________________________________
726 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
728 // Shape of the shower (see PHOS TDR)
729 // If you change this function, change also the gradient evaluation in ChiSquare()
731 //for the moment we neglect dependence on the incident angle.
733 Double_t r2 = x*x + z*z ;
734 Double_t r4 = r2*r2 ;
735 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
736 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
740 //____________________________________________________________________________
741 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
743 AliPHOSDigit ** maxAt,
744 Float_t * maxAtEnergy)
746 // Performs the unfolding of a cluster with nMax overlapping showers
748 Int_t nPar = 3 * nMax ;
749 Float_t * fitparameters = new Float_t[nPar] ;
751 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
754 // Fit failed, return and remove cluster
755 iniEmc->SetNExMax(-1) ;
756 delete[] fitparameters ;
760 // create ufolded rec points and fill them with new energy lists
761 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
762 // and later correct this number in acordance with actual energy deposition
764 Int_t nDigits = iniEmc->GetMultiplicity() ;
765 Float_t * efit = new Float_t[nDigits] ;
766 Float_t xDigit=0.,zDigit=0. ;
767 Float_t xpar=0.,zpar=0.,epar=0. ;
769 AliPHOSDigit * digit = 0 ;
770 Int_t * emcDigits = iniEmc->GetDigitsList() ;
776 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
777 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
778 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
779 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
783 while(iparam < nPar ){
784 xpar = fitparameters[iparam] ;
785 zpar = fitparameters[iparam+1] ;
786 epar = fitparameters[iparam+2] ;
788 // fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
789 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
790 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
794 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
795 // so that energy deposited in each cell is distributed betwin new clusters proportionally
796 // to its contribution to efit
798 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
802 while(iparam < nPar ){
803 xpar = fitparameters[iparam] ;
804 zpar = fitparameters[iparam+1] ;
805 epar = fitparameters[iparam+2] ;
807 // fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
809 AliPHOSEmcRecPoint * emcRP = 0 ;
811 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
813 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
814 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
816 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
817 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
818 fNumberOfEmcClusters++ ;
819 emcRP->SetNExMax((Int_t)nPar/3) ;
821 else{//create new entries in fCPVRecPoints
822 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
823 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
825 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
826 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
827 fNumberOfCpvClusters++ ;
831 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
832 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
833 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
834 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
835 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
836 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
837 eDigit = emcEnergies[iDigit] * ratio ;
838 emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId()) ) ;
842 delete[] fitparameters ;
847 //_____________________________________________________________________________
848 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
850 // Calculates the Chi square for the cluster unfolding minimization
851 // Number of parameters, Gradient, Chi squared, parameters, what to do
853 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
855 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
856 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
857 // A bit buggy way to get an access to the geometry
859 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
861 // TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
863 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
865 Int_t * emcDigits = emcRP->GetDigitsList() ;
867 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
869 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
876 for(iparam = 0 ; iparam < nPar ; iparam++)
877 Grad[iparam] = 0 ; // Will evaluate gradient
881 AliPHOSDigit * digit ;
884 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
886 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
892 geom->AbsToRelNumbering(digit->GetId(), relid) ;
894 geom->RelPosInModule(relid, xDigit, zDigit) ;
896 if(iflag == 2){ // calculate gradient
899 while(iParam < nPar ){
900 Double_t dx = (xDigit - x[iParam]) ;
902 Double_t dz = (zDigit - x[iParam]) ;
904 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
905 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
906 efit += x[iParam] * ShowerShape(dx,dz) ;
909 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
911 while(iParam < nPar ){
912 Double_t xpar = x[iParam] ;
913 Double_t zpar = x[iParam+1] ;
914 Double_t epar = x[iParam+2] ;
915 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
916 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
917 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
918 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
919 //DP: No incident angle dependence in gradient yet!!!!!!
920 Double_t r4 = dr*dr*dr*dr ;
921 Double_t r295 = TMath::Power(dr,2.95) ;
922 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
923 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
925 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
927 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
929 Grad[iParam] += shape ; // Derivative over energy
936 while(iparam < nPar ){
937 Double_t xpar = x[iparam] ;
938 Double_t zpar = x[iparam+1] ;
939 Double_t epar = x[iparam+2] ;
941 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
942 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
943 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
946 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
947 // Here we assume, that sigma = sqrt(E)
952 //____________________________________________________________________________
953 void AliPHOSClusterizerv1::Print(const Option_t *)const
955 // Print clusterizer parameters
958 TString taskName(GetName()) ;
959 taskName.ReplaceAll(Version(), "") ;
961 if( strcmp(GetName(), "") !=0 ) {
963 message = "\n--------------- %s %s -----------\n" ;
964 message += "Clusterizing digits from the file: %s\n" ;
965 message += " Branch: %s\n" ;
966 message += " EMC Clustering threshold = %f\n" ;
967 message += " EMC Local Maximum cut = %f\n" ;
968 message += " EMC Logarothmic weight = %f\n" ;
969 message += " CPV Clustering threshold = %f\n" ;
970 message += " CPV Local Maximum cut = %f\n" ;
971 message += " CPV Logarothmic weight = %f\n" ;
973 message += " Unfolding on\n" ;
975 message += " Unfolding off\n" ;
977 message += "------------------------------------------------------------------" ;
980 message = " AliPHOSClusterizerv1 not initialized " ;
982 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
987 fEmcClusteringThreshold,
990 fCpvClusteringThreshold,
994 //____________________________________________________________________________
995 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
997 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
999 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1000 fEMCRecPoints->GetEntriesFast(),
1001 fCPVRecPoints->GetEntriesFast() )) ;
1003 if(strstr(option,"all")) {
1004 printf("\n EMC clusters \n") ;
1005 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1007 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1008 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1010 rp->GetLocalPosition(locpos);
1012 rp->GetElipsAxis(lambda);
1015 primaries = rp->GetPrimaries(nprimaries);
1016 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1017 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1018 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1020 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1021 printf("%d ", primaries[iprimary] ) ;
1026 //Now plot CPV recPoints
1027 printf("\n CPV clusters \n") ;
1028 printf("Index Ene(MeV) Module X Y Z \n") ;
1029 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1030 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1033 rp->GetLocalPosition(locpos);
1035 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1036 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1037 locpos.X(), locpos.Y(), locpos.Z()) ;
1043 //____________________________________________________________________________
1044 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1046 //For each EMC rec. point set the distance to the nearest bad crystal.
1047 //Author: Boris Polichtchouk
1049 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1050 AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1053 fgCalibData->EmcBadChannelIds(badIds);
1055 AliPHOSEmcRecPoint* rp;
1058 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1059 TVector3 gposBadChannel; // global position of bad crystal
1062 Float_t dist,minDist;
1063 Int_t relid[4]={0,0,0,0} ;
1065 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1066 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1067 //evaluate distance to border
1068 relid[0]=rp->GetPHOSMod() ;
1071 Float_t xcorner,zcorner;
1072 fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1073 rp->GetLocalPosition(lpos) ;
1074 minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1075 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1076 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1077 if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1078 continue ; //bad channels can be in the module which does not exist in simulations.
1079 rp->GetGlobalPosition(gposRecPoint,gmat);
1080 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1081 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1082 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1083 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1084 dR = gposBadChannel-gposRecPoint;
1086 if(dist<minDist) minDist = dist;
1089 rp->SetDistanceToBadCrystal(minDist);
1093 //==================================================================================
1094 Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
1095 // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1097 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1099 //Determine rel.position of the cell absolute ID
1101 geom->AbsToRelNumbering(absId,relId);
1102 Int_t module=relId[0];
1103 Int_t row =relId[2];
1104 Int_t column=relId[3];
1106 Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1107 return amp*calibration ;
1110 Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1111 return amp*calibration ;
1114 //==================================================================================
1115 Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId)const{
1116 // Calibrate time in EMC digit
1118 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1120 //Determine rel.position of the cell absolute ID
1122 geom->AbsToRelNumbering(absId,relId);
1123 Int_t module=relId[0];
1124 Int_t row =relId[2];
1125 Int_t column=relId[3];
1130 time += fgCalibData->GetTimeShiftEmc(module,column,row);