1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
18 /* History of cvs commits:
21 * Revision 1.116 2007/10/01 20:24:08 kharlov
24 * Revision 1.115 2007/09/26 14:22:17 cvetan
25 * 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.
27 * Revision 1.114 2007/09/06 16:06:44 kharlov
28 * Absence of sorting results in loose of all unfolded clusters
30 * Revision 1.113 2007/08/28 12:55:07 policheh
31 * Loaders removed from the reconstruction code (C.Cheshkov)
33 * Revision 1.112 2007/08/22 09:20:50 hristov
34 * Updated QA classes (Yves)
36 * Revision 1.111 2007/08/08 12:11:28 kharlov
37 * Protection against uninitialized fQADM
39 * Revision 1.110 2007/08/07 14:16:00 kharlov
40 * Quality assurance added (Yves Schutz)
42 * Revision 1.109 2007/07/24 17:20:35 policheh
43 * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
44 * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
46 * Revision 1.108 2007/06/18 07:00:51 kharlov
47 * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
49 * Revision 1.107 2007/05/25 14:12:26 policheh
50 * Local to tracking CS transformation added for CPV rec. points
52 * Revision 1.106 2007/05/24 13:01:22 policheh
53 * Local to tracking CS transformation invoked for each EMC rec.point
55 * Revision 1.105 2007/05/02 13:41:22 kharlov
56 * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
58 * Revision 1.104 2007/04/27 16:55:53 kharlov
59 * Calibration stops if PHOS CDB objects do not exist
61 * Revision 1.103 2007/04/11 11:55:45 policheh
62 * SetDistancesToBadChannels() added.
64 * Revision 1.102 2007/03/28 19:18:15 kharlov
65 * RecPoints recalculation in TSM removed
67 * Revision 1.101 2007/03/06 06:51:27 kharlov
68 * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
70 * Revision 1.100 2007/01/10 11:07:26 kharlov
71 * Raw digits writing to file (B.Polichtchouk)
73 * Revision 1.99 2006/11/07 16:49:51 kharlov
74 * Corrections for next event switch in case of raw data (B.Polichtchouk)
76 * Revision 1.98 2006/10/27 17:14:27 kharlov
77 * Introduce AliDebug and AliLog (B.Polichtchouk)
79 * Revision 1.97 2006/08/29 11:41:19 kharlov
80 * Missing implementation of ctors and = operator are added
82 * Revision 1.96 2006/08/25 16:56:30 kharlov
83 * Compliance with Effective C++
85 * Revision 1.95 2006/08/11 12:36:26 cvetan
86 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
88 * Revision 1.94 2006/08/07 12:27:49 hristov
89 * Removing obsolete code which affected the event numbering scheme
91 * Revision 1.93 2006/08/01 12:20:17 cvetan
92 * 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
94 * Revision 1.92 2006/04/29 20:26:46 hristov
95 * Separate EMC and CPV calibration (Yu.Kharlov)
97 * Revision 1.91 2006/04/22 10:30:17 hristov
98 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
100 * Revision 1.90 2006/04/11 15:22:59 hristov
101 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
103 * Revision 1.89 2006/03/13 14:05:42 kharlov
104 * Calibration objects for EMC and CPV
106 * Revision 1.88 2006/01/11 08:54:52 hristov
107 * Additional protection in case no calibration entry was found
109 * Revision 1.87 2005/11/22 08:46:43 kharlov
110 * Updated to new CDB (Boris Polichtchouk)
112 * Revision 1.86 2005/11/14 21:52:43 hristov
115 * Revision 1.85 2005/09/27 16:08:08 hristov
116 * New version of CDB storage framework (A.Colla)
118 * Revision 1.84 2005/09/21 10:02:47 kharlov
119 * Reading calibration from CDB (Boris Polichtchouk)
121 * Revision 1.82 2005/09/02 15:43:13 kharlov
122 * Add comments in GetCalibrationParameters and Calibrate
124 * Revision 1.81 2005/09/02 14:32:07 kharlov
125 * Calibration of raw data
127 * Revision 1.80 2005/08/24 15:31:36 kharlov
128 * Setting raw digits flag
130 * Revision 1.79 2005/07/25 15:53:53 kharlov
133 * Revision 1.78 2005/05/28 14:19:04 schutz
134 * Compilation warnings fixed by T.P.
138 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
139 //////////////////////////////////////////////////////////////////////////////
140 // Clusterization class. Performs clusterization (collects neighbouring active cells) and
141 // unfolds the clusters having several local maxima.
142 // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
143 // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
144 // parameters including input digits branch title, thresholds etc.)
145 // This TTask is normally called from Reconstructioner, but can as well be used in
148 // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
149 // root [1] cl->Digits2Clusters(digitsTree,clusterTree)
150 // //finds RecPoints in the current event
151 // root [2] cl->SetDigitsBranch("digits2")
152 // //sets another title for Digitis (input) branch
153 // root [3] cl->SetRecPointsBranch("recp2")
154 // //sets another title four output branches
155 // root [4] cl->SetEmcLocalMaxCut(0.03)
156 // //set clusterization parameters
158 // --- ROOT system ---
163 #include "TBenchmark.h"
164 #include "TClonesArray.h"
166 // --- Standard library ---
168 // --- AliRoot header files ---
170 #include "AliConfig.h"
171 #include "AliPHOSGeometry.h"
172 #include "AliPHOSClusterizerv1.h"
173 #include "AliPHOSEmcRecPoint.h"
174 #include "AliPHOSCpvRecPoint.h"
175 #include "AliPHOSDigit.h"
176 #include "AliPHOSDigitizer.h"
177 #include "AliPHOSCalibrationDB.h"
178 #include "AliCDBManager.h"
179 #include "AliCDBStorage.h"
180 #include "AliCDBEntry.h"
181 #include "AliPHOSRecoParam.h"
182 #include "AliPHOSCalibData.h"
183 #include "AliPHOSReconstructor.h"
185 ClassImp(AliPHOSClusterizerv1)
187 //____________________________________________________________________________
188 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
189 AliPHOSClusterizer(),
190 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
191 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
192 fADCchanelEmc(0), fADCpedestalEmc(0),
193 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
194 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
195 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
196 fW0CPV(0), fEmcTimeGate(0),
199 // default ctor (to be used mainly by Streamer)
202 fDefaultInit = kTRUE ;
205 //____________________________________________________________________________
206 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
207 AliPHOSClusterizer(geom),
208 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
209 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
210 fADCchanelEmc(0), fADCpedestalEmc(0),
211 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
212 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
213 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
214 fW0CPV(0), fEmcTimeGate(0),
217 // ctor with the indication of the file where header Tree and digits Tree are stored
221 fDefaultInit = kFALSE ;
224 //____________________________________________________________________________
225 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
231 //____________________________________________________________________________
232 Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
234 // Convert EMC measured amplitude into real energy.
235 // Calibration parameters are taken from calibration data base for raw data,
236 // or from digitizer parameters for simulated data.
240 fGeom->AbsToRelNumbering(absId,relId) ;
241 Int_t module = relId[0];
242 Int_t column = relId[3];
243 Int_t row = relId[2];
244 if(absId <= fEmcCrystals) { // this is EMC
245 fADCchanelEmc = fgCalibData->GetADCchannelEmc (module,column,row);
246 return amp*fADCchanelEmc ;
250 if(absId <= fEmcCrystals) // this is EMC
251 return fADCpedestalEmc + amp*fADCchanelEmc ;
256 //____________________________________________________________________________
257 Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
259 // Convert digitized CPV amplitude into charge.
260 // Calibration parameters are taken from calibration data base for raw data,
261 // or from digitizer parameters for simulated data.
265 fGeom->AbsToRelNumbering(absId,relId) ;
266 Int_t module = relId[0];
267 Int_t column = relId[3];
268 Int_t row = relId[2];
269 if(absId > fEmcCrystals) { // this is CPV
270 fADCchanelCpv = fgCalibData->GetADCchannelCpv (module,column,row);
271 fADCpedestalCpv = fgCalibData->GetADCpedestalCpv(module,column,row);
272 return fADCpedestalCpv + amp*fADCchanelCpv ;
276 if(absId > fEmcCrystals) // this is CPV
277 return fADCpedestalCpv+ amp*fADCchanelCpv ;
282 //____________________________________________________________________________
283 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
285 // Steering method to perform clusterization for one event
286 // The input is the tree with digits.
287 // The output is the tree with clusters.
289 if(strstr(option,"tim"))
290 gBenchmark->Start("PHOSClusterizer");
292 if(strstr(option,"print")) {
299 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
300 fEMCRecPoints->GetEntries()));
301 if(AliLog::GetGlobalDebugLevel()>1)
302 fEMCRecPoints->Print();
309 if(strstr(option,"deb"))
310 PrintRecPoints(option) ;
312 if(strstr(option,"tim")){
313 gBenchmark->Stop("PHOSClusterizer");
314 AliInfo(Form("took %f seconds for Clusterizing\n",
315 gBenchmark->GetCpuTime("PHOSClusterizer")));
319 //____________________________________________________________________________
320 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
321 Int_t nPar, Float_t * fitparameters) const
323 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
324 // The initial values for fitting procedure are set equal to the positions of local maxima.
325 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
328 gMinuit->mncler(); // Reset Minuit's list of paramters
329 gMinuit->SetPrintLevel(-1) ; // No Printout
330 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
331 // To set the address of the minimization function
333 TList * toMinuit = new TList();
334 toMinuit->AddAt(emcRP,0) ;
335 toMinuit->AddAt(fDigitsArr,1) ;
336 toMinuit->AddAt(fGeom,2) ;
338 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
340 // filling initial values for fit parameters
341 AliPHOSDigit * digit ;
345 Int_t nDigits = (Int_t) nPar / 3 ;
349 for(iDigit = 0; iDigit < nDigits; iDigit++){
350 digit = maxAt[iDigit];
355 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
356 fGeom->RelPosInModule(relid, x, z) ;
358 Float_t energy = maxAtEnergy[iDigit] ;
360 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
363 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
366 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
369 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
372 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
375 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
380 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
385 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
386 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
387 gMinuit->SetMaxIterations(5);
388 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
390 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
392 if(ierflg == 4){ // Minimum not found
393 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
396 for(index = 0; index < nPar; index++){
399 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
400 fitparameters[index] = val ;
408 //____________________________________________________________________________
409 void AliPHOSClusterizerv1::GetCalibrationParameters()
411 // Set calibration parameters:
412 // if calibration database exists, they are read from database,
413 // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
415 // It is a user responsilibity to open CDB before reconstruction, for example:
416 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
419 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
420 if (fgCalibData->GetCalibDataEmc() == 0)
421 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
422 if (fgCalibData->GetCalibDataCpv() == 0)
423 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
427 //____________________________________________________________________________
428 void AliPHOSClusterizerv1::Init()
430 // Make all memory allocations which can not be done in default constructor.
431 // Attach the Clusterizer task to the list of PHOS tasks
433 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
436 gMinuit = new TMinuit(100);
438 GetCalibrationParameters() ;
442 //____________________________________________________________________________
443 void AliPHOSClusterizerv1::InitParameters()
446 fNumberOfCpvClusters = 0 ;
447 fNumberOfEmcClusters = 0 ;
449 const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
450 if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
452 const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv();
453 if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
455 fCpvClusteringThreshold = parCpv->GetClusteringThreshold();
456 fEmcClusteringThreshold = parEmc->GetClusteringThreshold();
458 fEmcLocMaxCut = parEmc->GetLocalMaxCut();
459 fCpvLocMaxCut = parCpv->GetLocalMaxCut();
461 fEmcMinE = parEmc->GetMinE();
462 fCpvMinE = parCpv->GetMinE();
464 fW0 = parEmc->GetLogWeight();
465 fW0CPV = parCpv->GetLogWeight();
467 fEmcTimeGate = 1.e-6 ;
473 fIsOldRCUFormat = kFALSE;
476 //____________________________________________________________________________
477 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
479 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
481 // = 2 are not neighbour but do not continue searching
482 // neighbours are defined as digits having at least a common vertex
483 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
484 // which is compared to a digit (d2) not yet in a cluster
489 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
492 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
494 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
495 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
496 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
498 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
499 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
500 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
504 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
505 rv = 2; // Difference in row numbers is too large to look further
511 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
518 //____________________________________________________________________________
519 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
521 // Remove digits with amplitudes below threshold.
522 // remove digits in bad channels
524 Bool_t isBadMap = 0 ;
525 if(fgCalibData->GetNumOfEmcBadChannels()){
530 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
531 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
532 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
533 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) ){
534 digits->RemoveAt(i) ;
537 if(isBadMap){ //check bad map now
539 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
540 if(fgCalibData->IsBadChannelEmc(relid[0],relid[2],relid[3])){
541 digits->RemoveAt(i) ;
547 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
548 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
549 digit->SetIndexInList(i) ;
553 //____________________________________________________________________________
554 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
556 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
560 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
562 if(digit->GetId() <= nEMC ) rv = kTRUE;
567 //____________________________________________________________________________
568 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
570 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
574 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
576 if(digit->GetId() > nEMC ) rv = kTRUE;
581 //____________________________________________________________________________
582 void AliPHOSClusterizerv1::WriteRecPoints()
585 // Creates new branches with given title
586 // fills and writes into TreeR.
589 //Evaluate position, dispersion and other RecPoint properties..
590 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
591 for(index = 0; index < nEmc; index++){
592 AliPHOSEmcRecPoint * rp =
593 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
594 rp->Purify(fEmcMinE) ;
595 if(rp->GetMultiplicity()==0){
596 fEMCRecPoints->RemoveAt(index) ;
601 // No vertex is available now, calculate corrections in PID
602 rp->EvalAll(fW0,fDigitsArr) ;
603 TVector3 fakeVtx(0.,0.,0.) ;
604 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
605 rp->EvalLocal2TrackingCSTransform();
607 fEMCRecPoints->Compress() ;
608 fEMCRecPoints->Sort() ;
609 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
610 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
611 dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
614 //For each rec.point set the distance to the nearest bad crystal (BVP)
615 SetDistancesToBadChannels();
617 //Now the same for CPV
618 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
619 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
620 rp->EvalAll(fW0CPV,fDigitsArr) ;
621 rp->EvalLocal2TrackingCSTransform();
623 fCPVRecPoints->Sort() ;
625 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
626 dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
628 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
630 if(fWrite){ //We write TreeR
635 //____________________________________________________________________________
636 void AliPHOSClusterizerv1::MakeClusters()
638 // Steering method to construct the clusters stored in a list of Reconstructed Points
639 // A cluster is defined as a list of neighbour digits
641 //Remove digits below threshold
642 CleanDigits(fDigitsArr) ;
644 TClonesArray * digitsC = static_cast<TClonesArray*>( fDigitsArr->Clone() ) ;
646 // Clusterization starts
648 TIter nextdigit(digitsC) ;
649 AliPHOSDigit * digit ;
650 Bool_t notremoved = kTRUE ;
652 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
655 AliPHOSRecPoint * clu = 0 ;
657 TArrayI clusterdigitslist(1500) ;
660 if (( IsInEmc (digit) &&
661 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
663 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
664 Int_t iDigitInCluster = 0 ;
666 if ( IsInEmc(digit) ) {
667 // start a new EMC RecPoint
668 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
669 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
671 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
672 clu = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
673 fNumberOfEmcClusters++ ;
674 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
675 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
677 digitsC->Remove(digit) ;
681 // start a new CPV cluster
682 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
683 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
685 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
687 clu = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
688 fNumberOfCpvClusters++ ;
689 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
690 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
692 digitsC->Remove(digit) ;
695 // Here we remove remaining EMC digits, which cannot make a cluster
698 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
700 digitsC->Remove(digit) ;
704 notremoved = kFALSE ;
711 AliPHOSDigit * digitN ;
713 while (index < iDigitInCluster){ // scan over digits already in cluster
714 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
716 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
717 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
719 case 0 : // not a neighbour
721 case 1 : // are neighbours
722 if (IsInEmc (digitN))
723 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
725 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
727 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
729 digitsC->Remove(digitN) ;
731 case 2 : // too far from each other
740 } // loop over cluster
751 //____________________________________________________________________________
752 void AliPHOSClusterizerv1::MakeUnfolding()
754 // Unfolds clusters using the shape of an ElectroMagnetic shower
755 // Performs unfolding of all EMC/CPV clusters
757 // Unfold first EMC clusters
758 if(fNumberOfEmcClusters > 0){
760 Int_t nModulesToUnfold = fGeom->GetNModules() ;
762 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
764 for(index = 0 ; index < numberofNotUnfolded ; index++){
766 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
767 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
770 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
771 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
772 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
773 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
775 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
776 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
777 fEMCRecPoints->Remove(emcRecPoint);
778 fEMCRecPoints->Compress() ;
780 fNumberOfEmcClusters -- ;
781 numberofNotUnfolded-- ;
784 emcRecPoint->SetNExMax(1) ; //Only one local maximum
788 delete[] maxAtEnergy ;
791 // Unfolding of EMC clusters finished
794 // Unfold now CPV clusters
795 if(fNumberOfCpvClusters > 0){
797 Int_t nModulesToUnfold = fGeom->GetNModules() ;
799 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
801 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
803 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
805 if(recPoint->GetPHOSMod()> nModulesToUnfold)
808 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
810 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
811 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
812 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
813 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
815 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
816 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
817 fCPVRecPoints->Remove(emcRecPoint);
818 fCPVRecPoints->Compress() ;
820 numberofCpvNotUnfolded-- ;
821 fNumberOfCpvClusters-- ;
825 delete[] maxAtEnergy ;
828 //Unfolding of Cpv clusters finished
832 //____________________________________________________________________________
833 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
835 // Shape of the shower (see PHOS TDR)
836 // If you change this function, change also the gradient evaluation in ChiSquare()
838 //for the moment we neglect dependence on the incident angle.
840 Double_t r2 = x*x + z*z ;
841 Double_t r4 = r2*r2 ;
842 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
843 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
847 //____________________________________________________________________________
848 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
850 AliPHOSDigit ** maxAt,
851 Float_t * maxAtEnergy)
853 // Performs the unfolding of a cluster with nMax overlapping showers
855 Int_t nPar = 3 * nMax ;
856 Float_t * fitparameters = new Float_t[nPar] ;
858 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
860 // Fit failed, return and remove cluster
861 iniEmc->SetNExMax(-1) ;
862 delete[] fitparameters ;
866 // create ufolded rec points and fill them with new energy lists
867 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
868 // and later correct this number in acordance with actual energy deposition
870 Int_t nDigits = iniEmc->GetMultiplicity() ;
871 Float_t * efit = new Float_t[nDigits] ;
872 Float_t xDigit=0.,zDigit=0. ;
873 Float_t xpar=0.,zpar=0.,epar=0. ;
875 AliPHOSDigit * digit = 0 ;
876 Int_t * emcDigits = iniEmc->GetDigitsList() ;
882 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
883 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
884 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
885 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
889 while(iparam < nPar ){
890 xpar = fitparameters[iparam] ;
891 zpar = fitparameters[iparam+1] ;
892 epar = fitparameters[iparam+2] ;
894 // fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
895 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
896 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
901 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
902 // so that energy deposited in each cell is distributed betwin new clusters proportionally
903 // to its contribution to efit
905 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
909 while(iparam < nPar ){
910 xpar = fitparameters[iparam] ;
911 zpar = fitparameters[iparam+1] ;
912 epar = fitparameters[iparam+2] ;
914 // fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
916 AliPHOSEmcRecPoint * emcRP = 0 ;
918 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
920 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
921 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
923 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
924 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
925 fNumberOfEmcClusters++ ;
926 emcRP->SetNExMax((Int_t)nPar/3) ;
928 else{//create new entries in fCPVRecPoints
929 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
930 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
932 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
933 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
934 fNumberOfCpvClusters++ ;
938 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
939 digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
940 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
941 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
942 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
943 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
944 eDigit = emcEnergies[iDigit] * ratio ;
945 emcRP->AddDigit( *digit, eDigit ) ;
949 delete[] fitparameters ;
954 //_____________________________________________________________________________
955 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
957 // Calculates the Chi square for the cluster unfolding minimization
958 // Number of parameters, Gradient, Chi squared, parameters, what to do
960 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
962 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
963 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
964 // A bit buggy way to get an access to the geometry
966 AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
968 // TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
970 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
972 Int_t * emcDigits = emcRP->GetDigitsList() ;
974 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
976 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
983 for(iparam = 0 ; iparam < nPar ; iparam++)
984 Grad[iparam] = 0 ; // Will evaluate gradient
988 AliPHOSDigit * digit ;
991 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
993 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
999 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1001 geom->RelPosInModule(relid, xDigit, zDigit) ;
1003 if(iflag == 2){ // calculate gradient
1006 while(iParam < nPar ){
1007 Double_t dx = (xDigit - x[iParam]) ;
1009 Double_t dz = (zDigit - x[iParam]) ;
1011 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
1012 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
1013 efit += x[iParam] * ShowerShape(dx,dz) ;
1016 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
1018 while(iParam < nPar ){
1019 Double_t xpar = x[iParam] ;
1020 Double_t zpar = x[iParam+1] ;
1021 Double_t epar = x[iParam+2] ;
1022 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1023 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1024 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1025 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1026 //DP: No incident angle dependence in gradient yet!!!!!!
1027 Double_t r4 = dr*dr*dr*dr ;
1028 Double_t r295 = TMath::Power(dr,2.95) ;
1029 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1030 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1032 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1034 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1036 Grad[iParam] += shape ; // Derivative over energy
1043 while(iparam < nPar ){
1044 Double_t xpar = x[iparam] ;
1045 Double_t zpar = x[iparam+1] ;
1046 Double_t epar = x[iparam+2] ;
1048 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
1049 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
1050 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
1053 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1054 // Here we assume, that sigma = sqrt(E)
1059 //____________________________________________________________________________
1060 void AliPHOSClusterizerv1::Print(const Option_t *)const
1062 // Print clusterizer parameters
1065 TString taskName(GetName()) ;
1066 taskName.ReplaceAll(Version(), "") ;
1068 if( strcmp(GetName(), "") !=0 ) {
1070 message = "\n--------------- %s %s -----------\n" ;
1071 message += "Clusterizing digits from the file: %s\n" ;
1072 message += " Branch: %s\n" ;
1073 message += " EMC Clustering threshold = %f\n" ;
1074 message += " EMC Local Maximum cut = %f\n" ;
1075 message += " EMC Logarothmic weight = %f\n" ;
1076 message += " CPV Clustering threshold = %f\n" ;
1077 message += " CPV Local Maximum cut = %f\n" ;
1078 message += " CPV Logarothmic weight = %f\n" ;
1080 message += " Unfolding on\n" ;
1082 message += " Unfolding off\n" ;
1084 message += "------------------------------------------------------------------" ;
1087 message = " AliPHOSClusterizerv1 not initialized " ;
1089 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1094 fEmcClusteringThreshold,
1097 fCpvClusteringThreshold,
1101 //____________________________________________________________________________
1102 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1104 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1106 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1107 fEMCRecPoints->GetEntriesFast(),
1108 fCPVRecPoints->GetEntriesFast() )) ;
1110 if(strstr(option,"all")) {
1111 printf("\n EMC clusters \n") ;
1112 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1114 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1115 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1117 rp->GetLocalPosition(locpos);
1119 rp->GetElipsAxis(lambda);
1122 primaries = rp->GetPrimaries(nprimaries);
1123 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1124 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1125 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1127 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1128 printf("%d ", primaries[iprimary] ) ;
1133 //Now plot CPV recPoints
1134 printf("\n CPV clusters \n") ;
1135 printf("Index Ene(MeV) Module X Y Z \n") ;
1136 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1137 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1140 rp->GetLocalPosition(locpos);
1142 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1143 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1144 locpos.X(), locpos.Y(), locpos.Z()) ;
1150 //____________________________________________________________________________
1151 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1153 //For each EMC rec. point set the distance to the nearest bad crystal.
1154 //Author: Boris Polichtchouk
1156 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1157 AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels()));
1160 fgCalibData->EmcBadChannelIds(badIds);
1162 AliPHOSEmcRecPoint* rp;
1165 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1166 TVector3 gposBadChannel; // global position of bad crystal
1169 Float_t dist,minDist;
1171 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1172 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1175 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1176 rp->GetGlobalPosition(gposRecPoint,gmat);
1177 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1178 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1179 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1180 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1181 dR = gposBadChannel-gposRecPoint;
1183 if(dist<minDist) minDist = dist;
1186 rp->SetDistanceToBadCrystal(minDist);