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),
201 fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.),
204 // default ctor (to be used mainly by Streamer)
206 fDefaultInit = kTRUE ;
208 for(Int_t i=0; i<53760; i++){
213 //____________________________________________________________________________
214 AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
215 AliPHOSClusterizer(geom),
216 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
218 fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
219 fEmcClusteringThreshold(0), fCpvClusteringThreshold(0),
220 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
222 fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.),
225 // ctor with the indication of the file where header Tree and digits Tree are stored
227 for(Int_t i=0; i<53760; i++){
232 fDefaultInit = kFALSE ;
235 //____________________________________________________________________________
236 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
241 //____________________________________________________________________________
242 void AliPHOSClusterizerv1::Digits2Clusters(Option_t *option)
244 // Steering method to perform clusterization for one event
245 // The input is the tree with digits.
246 // The output is the tree with clusters.
248 if(strstr(option,"tim"))
249 gBenchmark->Start("PHOSClusterizer");
251 if(strstr(option,"print")) {
258 AliDebug(2,Form(" ---- Printing clusters (%d)\n",
259 fEMCRecPoints->GetEntries()));
260 if(AliLog::GetGlobalDebugLevel()>1)
261 fEMCRecPoints->Print();
268 if(strstr(option,"deb"))
269 PrintRecPoints(option) ;
271 if(strstr(option,"tim")){
272 gBenchmark->Stop("PHOSClusterizer");
273 AliInfo(Form("took %f seconds for Clusterizing\n",
274 gBenchmark->GetCpuTime("PHOSClusterizer")));
276 fEMCRecPoints->Delete();
277 fCPVRecPoints->Delete();
280 //____________________________________________________________________________
281 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
282 Int_t nPar, Float_t * fitparameters) const
284 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
285 // The initial values for fitting procedure are set equal to the positions of local maxima.
286 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
289 if(!gMinuit) //it was deleted by someone else
290 gMinuit = new TMinuit(100) ;
291 gMinuit->mncler(); // Reset Minuit's list of paramters
292 gMinuit->SetPrintLevel(-1) ; // No Printout
293 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
294 // To set the address of the minimization function
296 TList * toMinuit = new TList();
297 toMinuit->AddAt(emcRP,0) ;
298 toMinuit->AddAt(fDigitsArr,1) ;
299 toMinuit->AddAt(fGeom,2) ;
301 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
303 // filling initial values for fit parameters
304 AliPHOSDigit * digit ;
308 Int_t nDigits = (Int_t) nPar / 3 ;
312 for(iDigit = 0; iDigit < nDigits; iDigit++){
313 digit = maxAt[iDigit];
318 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
319 fGeom->RelPosInModule(relid, x, z) ;
321 Float_t energy = maxAtEnergy[iDigit] ;
323 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
326 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
329 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
332 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
335 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
338 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
343 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
348 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
349 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
350 gMinuit->SetMaxIterations(5);
351 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
353 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
355 if(ierflg == 4){ // Minimum not found
356 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
359 for(index = 0; index < nPar; index++){
362 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
363 fitparameters[index] = val ;
372 //____________________________________________________________________________
373 void AliPHOSClusterizerv1::Init()
375 // Make all memory allocations which can not be done in default constructor.
376 // Attach the Clusterizer task to the list of PHOS tasks
378 fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
381 gMinuit = new TMinuit(100);
384 fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
385 if (fgCalibData->GetCalibDataEmc() == 0)
386 AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
387 if (fgCalibData->GetCalibDataCpv() == 0)
388 AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
392 //____________________________________________________________________________
393 void AliPHOSClusterizerv1::InitParameters()
396 fNumberOfCpvClusters = 0 ;
397 fNumberOfEmcClusters = 0 ;
399 const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam();
400 if(!recoParam) AliFatal("Reconstruction parameters are not set!");
404 fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold();
405 fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold();
407 fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut();
408 fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut();
410 fW0 = recoParam->GetEMCLogWeight();
411 fW0CPV = recoParam->GetCPVLogWeight();
413 fTimeGateLowAmp = recoParam->GetTimeGateAmpThresh() ;
414 fTimeGateLow = recoParam->GetTimeGateLow() ;
415 fTimeGateHigh = recoParam->GetTimeGateHigh() ;
417 fEcoreRadius = recoParam->GetEMCEcoreRadius();
419 fToUnfold = recoParam->EMCToUnfold() ;
424 //____________________________________________________________________________
425 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
427 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
429 // = 2 are not neighbour but do not continue searching
430 // =-1 are not neighbour, continue searching, but do not look before d2 next time
431 // neighbours are defined as digits having at least a common vertex
432 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
433 // which is compared to a digit (d2) not yet in a cluster
436 fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
439 fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
441 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
442 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
443 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
445 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex
446 // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side
447 if((relid1[1] != 0) || CheckTimeGate(d1->GetTime(),d1->GetEnergy(),d2->GetTime(),d2->GetEnergy()))
451 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
452 return 2; // Difference in row numbers is too large to look further
458 if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module
460 if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1)
469 //____________________________________________________________________________
470 Bool_t AliPHOSClusterizerv1::CheckTimeGate(Float_t t1, Float_t amp1, Float_t t2, Float_t amp2)const{
471 //Check if two cells have reasonable time difference
472 //Note that at low amplitude time is defined up to 1 tick == 100 ns.
473 if(amp1<fTimeGateLowAmp || amp2<fTimeGateLowAmp){
474 return (TMath::Abs(t1 - t2 ) < fTimeGateLow) ;
476 else{ //Time should be measured with good accuracy
477 return (TMath::Abs(t1 - t2 ) < fTimeGateHigh) ;
481 //____________________________________________________________________________
482 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
484 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
488 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
490 if(digit->GetId() <= nEMC ) rv = kTRUE;
495 //____________________________________________________________________________
496 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
498 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
502 Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
504 if(digit->GetId() > nEMC ) rv = kTRUE;
509 //____________________________________________________________________________
510 void AliPHOSClusterizerv1::WriteRecPoints()
513 // Creates new branches with given title
514 // fills and writes into TreeR.
517 //Evaluate position, dispersion and other RecPoint properties..
518 Int_t nEmc = fEMCRecPoints->GetEntriesFast();
519 Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy
520 TVector3 fakeVtx(0.,0.,0.) ;
521 for(index = 0; index < nEmc; index++){
522 AliPHOSEmcRecPoint * rp =
523 static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
524 rp->Purify(emcMinE,fDigitsArr) ;
525 if(rp->GetMultiplicity()==0){
526 fEMCRecPoints->RemoveAt(index) ;
531 // No vertex is available now, calculate corrections in PID
532 rp->EvalAll(fDigitsArr) ;
533 rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ;
534 rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
535 rp->EvalLocal2TrackingCSTransform();
537 fEMCRecPoints->Compress() ;
538 fEMCRecPoints->Sort() ;
539 // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
540 for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
541 static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
544 //For each rec.point set the distance to the nearest bad crystal (BVP)
545 SetDistancesToBadChannels();
547 //Now the same for CPV
548 for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
549 AliPHOSCpvRecPoint * rp = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
550 rp->EvalAll(fDigitsArr) ;
551 rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ;
552 rp->EvalLocal2TrackingCSTransform();
554 fCPVRecPoints->Sort() ;
556 for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
557 static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
559 fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
561 if(fWrite){ //We write TreeR
566 //____________________________________________________________________________
567 void AliPHOSClusterizerv1::MakeClusters()
569 // Steering method to construct the clusters stored in a list of Reconstructed Points
570 // A cluster is defined as a list of neighbour digits
572 fNumberOfCpvClusters = 0 ;
573 fNumberOfEmcClusters = 0 ;
575 //Mark all digits as unused yet
576 const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;)
577 Int_t nDigits=fDigitsArr->GetEntriesFast() ;
579 for(Int_t i=0; i<nDigits; i++){
582 Int_t iFirst = 0 ; //first index of digit which potentially can be a part of cluster
583 //e.g. first digit in this module, first CPV digit etc.
584 AliPHOSDigit * digit ;
585 TArrayI clusterdigitslist(maxNDigits) ;
586 AliPHOSRecPoint * clu = 0 ;
587 for(Int_t i=0; i<nDigits; i++){
591 digit=static_cast<AliPHOSDigit*>(fDigitsArr->At(i)) ;
597 //is this digit so energetic that start cluster?
598 if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
599 ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) {
600 Int_t iDigitInCluster = 0 ;
601 if ( IsInEmc(digit) ) {
602 // start a new EMC RecPoint
603 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
604 fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
606 fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
607 clu = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
608 fNumberOfEmcClusters++ ;
609 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG())) ;
610 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
612 fDigitsUsed[i]=kTRUE ;
614 // start a new CPV cluster
615 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
616 fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
618 fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
619 clu = static_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
620 fNumberOfCpvClusters++ ;
621 clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV
622 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
624 fDigitsUsed[i]=kTRUE ;
627 //Now scan remaining digits in list to find neigbours of our seed
629 AliPHOSDigit * digitN ;
631 while (index < iDigitInCluster){ // scan over digits already in cluster
632 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(clusterdigitslist[index]) ) ;
634 for(Int_t j=iFirst; j<nDigits; j++){
635 if (iDigitInCluster >= maxNDigits) {
636 AliError(Form("The number of digits in cluster is more than %d, skip the rest of event",
641 continue ; //look through remaining digits
642 digitN = static_cast<AliPHOSDigit*>( fDigitsArr->At(j) ) ;
643 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
645 case -1: //too early (e.g. previous module), do not look before j at subsequent passes
648 case 0 : // not a neighbour
650 case 1 : // are neighbours
651 clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digitN->GetId()),CalibrateT(digitN->GetTime(),digitN->GetId(),digit->IsLG())) ;
652 clusterdigitslist[iDigitInCluster] = j ;
654 fDigitsUsed[j]=kTRUE ;
656 case 2 : // too far from each other
662 endOfLoop: ; //scanned all possible neighbours for this digit
664 } // loop over cluster
670 //____________________________________________________________________________
671 void AliPHOSClusterizerv1::MakeUnfolding()
673 // Unfolds clusters using the shape of an ElectroMagnetic shower
674 // Performs unfolding of all EMC/CPV clusters
676 // Unfold first EMC clusters
677 if(fNumberOfEmcClusters > 0){
679 Int_t nModulesToUnfold = fGeom->GetNModules() ;
681 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
683 for(index = 0 ; index < numberofNotUnfolded ; index++){
685 AliPHOSEmcRecPoint * emcRecPoint = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
686 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
689 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
690 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
691 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
692 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
694 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
695 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
697 fEMCRecPoints->Remove(emcRecPoint);
698 fEMCRecPoints->Compress() ;
700 fNumberOfEmcClusters -- ;
701 numberofNotUnfolded-- ;
704 emcRecPoint->SetNExMax(1) ; //Only one local maximum
708 delete[] maxAtEnergy ;
711 // Unfolding of EMC clusters finished
714 // Unfold now CPV clusters
715 if(fNumberOfCpvClusters > 0){
717 Int_t nModulesToUnfold = fGeom->GetNModules() ;
719 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
721 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
723 AliPHOSRecPoint * recPoint = static_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
725 if(recPoint->GetPHOSMod()> nModulesToUnfold)
728 AliPHOSEmcRecPoint * emcRecPoint = static_cast<AliPHOSEmcRecPoint*>(recPoint) ;
730 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
731 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
732 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
733 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
735 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
736 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
737 fCPVRecPoints->Remove(emcRecPoint);
738 fCPVRecPoints->Compress() ;
740 numberofCpvNotUnfolded-- ;
741 fNumberOfCpvClusters-- ;
745 delete[] maxAtEnergy ;
748 //Unfolding of Cpv clusters finished
752 //____________________________________________________________________________
753 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
755 // Shape of the shower (see PHOS TDR)
756 // If you change this function, change also the gradient evaluation in ChiSquare()
758 //for the moment we neglect dependence on the incident angle.
760 Double_t r2 = x*x + z*z ;
761 Double_t r4 = r2*r2 ;
762 Double_t r295 = TMath::Power(r2, 2.95/2.) ;
763 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
767 //____________________________________________________________________________
768 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
770 AliPHOSDigit ** maxAt,
771 Float_t * maxAtEnergy)
773 // Performs the unfolding of a cluster with nMax overlapping showers
775 Int_t nPar = 3 * nMax ;
776 Float_t * fitparameters = new Float_t[nPar] ;
778 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
781 // Fit failed, return and remove cluster
782 iniEmc->SetNExMax(-1) ;
783 delete[] fitparameters ;
787 // create ufolded rec points and fill them with new energy lists
788 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
789 // and later correct this number in acordance with actual energy deposition
791 Int_t nDigits = iniEmc->GetMultiplicity() ;
792 Float_t * efit = new Float_t[nDigits] ;
793 Float_t xDigit=0.,zDigit=0. ;
794 Float_t xpar=0.,zpar=0.,epar=0. ;
796 AliPHOSDigit * digit = 0 ;
797 Int_t * emcDigits = iniEmc->GetDigitsList() ;
803 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
804 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
805 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
806 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
810 while(iparam < nPar ){
811 xpar = fitparameters[iparam] ;
812 zpar = fitparameters[iparam+1] ;
813 epar = fitparameters[iparam+2] ;
815 // fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
816 // efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
817 efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
821 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
822 // so that energy deposited in each cell is distributed betwin new clusters proportionally
823 // to its contribution to efit
825 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
829 while(iparam < nPar ){
830 xpar = fitparameters[iparam] ;
831 zpar = fitparameters[iparam+1] ;
832 epar = fitparameters[iparam+2] ;
834 // fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
836 AliPHOSEmcRecPoint * emcRP = 0 ;
838 if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
840 if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
841 fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
843 (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
844 emcRP = static_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
845 fNumberOfEmcClusters++ ;
846 emcRP->SetNExMax((Int_t)nPar/3) ;
848 else{//create new entries in fCPVRecPoints
849 if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
850 fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
852 (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
853 emcRP = static_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
854 fNumberOfCpvClusters++ ;
858 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
859 digit = static_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
860 fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
861 fGeom->RelPosInModule(relid, xDigit, zDigit) ;
862 // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
863 ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
864 eDigit = emcEnergies[iDigit] * ratio ;
865 emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG()) ) ;
869 delete[] fitparameters ;
874 //_____________________________________________________________________________
875 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
877 // Calculates the Chi square for the cluster unfolding minimization
878 // Number of parameters, Gradient, Chi squared, parameters, what to do
880 TList * toMinuit = static_cast<TList*>( gMinuit->GetObjectFit() ) ;
882 AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
883 TClonesArray * digits = static_cast<TClonesArray*>( toMinuit->At(1) ) ;
884 // A bit buggy way to get an access to the geometry
886 AliPHOSGeometry *geom = static_cast<AliPHOSGeometry *>(toMinuit->At(2));
888 // TVector3 * vtx = static_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
890 // AliPHOSEmcRecPoint * emcRP = static_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
892 Int_t * emcDigits = emcRP->GetDigitsList() ;
894 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
896 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
903 for(iparam = 0 ; iparam < nPar ; iparam++)
904 Grad[iparam] = 0 ; // Will evaluate gradient
908 AliPHOSDigit * digit ;
911 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
913 digit = static_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
919 geom->AbsToRelNumbering(digit->GetId(), relid) ;
921 geom->RelPosInModule(relid, xDigit, zDigit) ;
923 if(iflag == 2){ // calculate gradient
926 while(iParam < nPar ){
927 Double_t dx = (xDigit - x[iParam]) ;
929 Double_t dz = (zDigit - x[iParam]) ;
931 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
932 // efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
933 efit += x[iParam] * ShowerShape(dx,dz) ;
936 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
938 while(iParam < nPar ){
939 Double_t xpar = x[iParam] ;
940 Double_t zpar = x[iParam+1] ;
941 Double_t epar = x[iParam+2] ;
942 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
943 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
944 // Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
945 Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
946 //DP: No incident angle dependence in gradient yet!!!!!!
947 Double_t r4 = dr*dr*dr*dr ;
948 Double_t r295 = TMath::Power(dr,2.95) ;
949 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
950 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
952 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
954 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
956 Grad[iParam] += shape ; // Derivative over energy
963 while(iparam < nPar ){
964 Double_t xpar = x[iparam] ;
965 Double_t zpar = x[iparam+1] ;
966 Double_t epar = x[iparam+2] ;
968 // fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
969 // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
970 efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
973 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
974 // Here we assume, that sigma = sqrt(E)
979 //____________________________________________________________________________
980 void AliPHOSClusterizerv1::Print(const Option_t *)const
982 // Print clusterizer parameters
985 TString taskName(GetName()) ;
986 taskName.ReplaceAll(Version(), "") ;
988 if( strcmp(GetName(), "") !=0 ) {
990 message = "\n--------------- %s %s -----------\n" ;
991 message += "Clusterizing digits from the file: %s\n" ;
992 message += " Branch: %s\n" ;
993 message += " EMC Clustering threshold = %f\n" ;
994 message += " EMC Local Maximum cut = %f\n" ;
995 message += " EMC Logarothmic weight = %f\n" ;
996 message += " CPV Clustering threshold = %f\n" ;
997 message += " CPV Local Maximum cut = %f\n" ;
998 message += " CPV Logarothmic weight = %f\n" ;
1000 message += " Unfolding on\n" ;
1002 message += " Unfolding off\n" ;
1004 message += "------------------------------------------------------------------" ;
1007 message = " AliPHOSClusterizerv1 not initialized " ;
1009 AliInfo(Form("%s, %s %s %s %s %f %f %f %f %f %f", message.Data(),
1014 fEmcClusteringThreshold,
1017 fCpvClusteringThreshold,
1021 //____________________________________________________________________________
1022 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1024 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1026 AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
1027 fEMCRecPoints->GetEntriesFast(),
1028 fCPVRecPoints->GetEntriesFast() )) ;
1030 if(strstr(option,"all")) {
1031 printf("\n EMC clusters \n") ;
1032 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1034 for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
1035 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
1037 rp->GetLocalPosition(locpos);
1039 rp->GetElipsAxis(lambda);
1042 primaries = rp->GetPrimaries(nprimaries);
1043 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1044 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1045 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1047 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1048 printf("%d ", primaries[iprimary] ) ;
1053 //Now plot CPV recPoints
1054 printf("\n CPV clusters \n") ;
1055 printf("Index Ene(MeV) Module X Y Z \n") ;
1056 for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
1057 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
1060 rp->GetLocalPosition(locpos);
1062 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1063 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1064 locpos.X(), locpos.Y(), locpos.Z()) ;
1070 //____________________________________________________________________________
1071 void AliPHOSClusterizerv1::SetDistancesToBadChannels()
1073 //For each EMC rec. point set the distance to the nearest bad crystal.
1074 //Author: Boris Polichtchouk
1076 if(!fgCalibData->GetNumOfEmcBadChannels()) return;
1079 memset(badIds,0,8000*sizeof(Int_t));
1080 fgCalibData->EmcBadChannelIds(badIds);
1082 AliPHOSEmcRecPoint* rp;
1085 TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
1086 TVector3 gposBadChannel; // global position of bad crystal
1089 Float_t dist,minDist;
1090 Int_t relid[4]={0,0,0,0} ;
1092 for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
1093 rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
1094 //evaluate distance to border
1095 relid[0]=rp->GetPHOSMod() ;
1098 Float_t xcorner,zcorner;
1099 fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell
1100 rp->GetLocalPosition(lpos) ;
1101 minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size
1102 for(Int_t iBad=0; iBad<fgCalibData->GetNumOfEmcBadChannels(); iBad++) {
1103 fGeom->AbsToRelNumbering(badIds[iBad],relid) ;
1104 if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since
1105 continue ; //bad channels can be in the module which does not exist in simulations.
1106 rp->GetGlobalPosition(gposRecPoint,gmat);
1107 fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
1108 AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
1109 gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
1110 gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
1111 dR = gposBadChannel-gposRecPoint;
1113 if(dist<minDist) minDist = dist;
1116 rp->SetDistanceToBadCrystal(minDist);
1120 //==================================================================================
1121 Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{
1122 // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB
1124 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1126 //Determine rel.position of the cell absolute ID
1128 geom->AbsToRelNumbering(absId,relId);
1129 Int_t module=relId[0];
1130 Int_t row =relId[2];
1131 Int_t column=relId[3];
1133 Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row);
1134 return amp*calibration ;
1137 Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row);
1138 return amp*calibration ;
1141 //==================================================================================
1142 Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId,Bool_t isLG)const{
1143 // Calibrate time in EMC digit
1145 const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ;
1147 //Determine rel.position of the cell absolute ID
1149 geom->AbsToRelNumbering(absId,relId);
1150 Int_t module=relId[0];
1151 Int_t row =relId[2];
1152 Int_t column=relId[3];
1158 time += fgCalibData->GetLGTimeShiftEmc(module,column,row);
1160 time += fgCalibData->GetTimeShiftEmc(module,column,row);