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.99 2006/11/07 16:49:51 kharlov
22 * Corrections for next event switch in case of raw data (B.Polichtchouk)
24 * Revision 1.98 2006/10/27 17:14:27 kharlov
25 * Introduce AliDebug and AliLog (B.Polichtchouk)
27 * Revision 1.97 2006/08/29 11:41:19 kharlov
28 * Missing implementation of ctors and = operator are added
30 * Revision 1.96 2006/08/25 16:56:30 kharlov
31 * Compliance with Effective C++
33 * Revision 1.95 2006/08/11 12:36:26 cvetan
34 * Update of the PHOS code needed in order to read and reconstruct the beam test raw data (i.e. without an existing galice.root)
36 * Revision 1.94 2006/08/07 12:27:49 hristov
37 * Removing obsolete code which affected the event numbering scheme
39 * Revision 1.93 2006/08/01 12:20:17 cvetan
40 * 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
42 * Revision 1.92 2006/04/29 20:26:46 hristov
43 * Separate EMC and CPV calibration (Yu.Kharlov)
45 * Revision 1.91 2006/04/22 10:30:17 hristov
46 * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
48 * Revision 1.90 2006/04/11 15:22:59 hristov
49 * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
51 * Revision 1.89 2006/03/13 14:05:42 kharlov
52 * Calibration objects for EMC and CPV
54 * Revision 1.88 2006/01/11 08:54:52 hristov
55 * Additional protection in case no calibration entry was found
57 * Revision 1.87 2005/11/22 08:46:43 kharlov
58 * Updated to new CDB (Boris Polichtchouk)
60 * Revision 1.86 2005/11/14 21:52:43 hristov
63 * Revision 1.85 2005/09/27 16:08:08 hristov
64 * New version of CDB storage framework (A.Colla)
66 * Revision 1.84 2005/09/21 10:02:47 kharlov
67 * Reading calibration from CDB (Boris Polichtchouk)
69 * Revision 1.82 2005/09/02 15:43:13 kharlov
70 * Add comments in GetCalibrationParameters and Calibrate
72 * Revision 1.81 2005/09/02 14:32:07 kharlov
73 * Calibration of raw data
75 * Revision 1.80 2005/08/24 15:31:36 kharlov
76 * Setting raw digits flag
78 * Revision 1.79 2005/07/25 15:53:53 kharlov
81 * Revision 1.78 2005/05/28 14:19:04 schutz
82 * Compilation warnings fixed by T.P.
86 //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
87 //////////////////////////////////////////////////////////////////////////////
88 // Clusterization class. Performs clusterization (collects neighbouring active cells) and
89 // unfolds the clusters having several local maxima.
90 // Results are stored in TreeR#, branches PHOSEmcRP (EMC recPoints),
91 // PHOSCpvRP (CPV RecPoints) and AliPHOSClusterizer (Clusterizer with all
92 // parameters including input digits branch title, thresholds etc.)
93 // This TTask is normally called from Reconstructioner, but can as well be used in
96 // root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
97 // Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
98 // // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
99 // // and saves recpoints in branch named "recpointsname" (default = "digitsname")
100 // root [1] cl->ExecuteTask()
101 // //finds RecPoints in all events stored in galice.root
102 // root [2] cl->SetDigitsBranch("digits2")
103 // //sets another title for Digitis (input) branch
104 // root [3] cl->SetRecPointsBranch("recp2")
105 // //sets another title four output branches
106 // root [4] cl->SetEmcLocalMaxCut(0.03)
107 // //set clusterization parameters
108 // root [5] cl->ExecuteTask("deb all time")
109 // //once more finds RecPoints options are
110 // // deb - print number of found rec points
111 // // deb all - print number of found RecPoints and some their characteristics
112 // // time - print benchmarking results
114 // --- ROOT system ---
119 #include "TBenchmark.h"
121 // --- Standard library ---
123 // --- AliRoot header files ---
125 #include "AliRunLoader.h"
126 #include "AliPHOSGetter.h"
127 #include "AliPHOSGeometry.h"
128 #include "AliPHOSClusterizerv1.h"
129 #include "AliPHOSEmcRecPoint.h"
130 #include "AliPHOSCpvRecPoint.h"
131 #include "AliPHOSDigit.h"
132 #include "AliPHOSDigitizer.h"
133 #include "AliPHOSCalibrationDB.h"
134 #include "AliCDBManager.h"
135 #include "AliCDBStorage.h"
136 #include "AliCDBEntry.h"
138 ClassImp(AliPHOSClusterizerv1)
140 //____________________________________________________________________________
141 AliPHOSClusterizerv1::AliPHOSClusterizerv1() :
142 AliPHOSClusterizer(),
143 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
144 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
145 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
146 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
147 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
148 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
149 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
152 // default ctor (to be used mainly by Streamer)
155 fDefaultInit = kTRUE ;
158 //____________________________________________________________________________
159 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
160 AliPHOSClusterizer(alirunFileName, eventFolderName),
161 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
162 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
163 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
164 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
165 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
166 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
167 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
170 // ctor with the indication of the file where header Tree and digits Tree are stored
174 fDefaultInit = kFALSE ;
177 //____________________________________________________________________________
178 AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
179 AliPHOSClusterizer(obj),
180 fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
181 fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
182 fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
183 fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
184 fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
185 fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
186 fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
191 //____________________________________________________________________________
192 AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
197 //____________________________________________________________________________
198 const TString AliPHOSClusterizerv1::BranchName() const
203 //____________________________________________________________________________
204 Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
206 // Convert EMC measured amplitude into real energy.
207 // Calibration parameters are taken from calibration data base for raw data,
208 // or from digitizer parameters for simulated data.
212 AliPHOSGetter *gime = AliPHOSGetter::Instance();
213 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
214 Int_t module = relId[0];
215 Int_t column = relId[3];
216 Int_t row = relId[2];
217 if(absId <= fEmcCrystals) { // this is EMC
218 fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
219 return amp*fADCchanelEmc ;
223 if(absId <= fEmcCrystals) // this is EMC
224 return fADCpedestalEmc + amp*fADCchanelEmc ;
229 //____________________________________________________________________________
230 Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
232 // Convert digitized CPV amplitude into charge.
233 // Calibration parameters are taken from calibration data base for raw data,
234 // or from digitizer parameters for simulated data.
238 AliPHOSGetter *gime = AliPHOSGetter::Instance();
239 gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
240 Int_t module = relId[0];
241 Int_t column = relId[3];
242 Int_t row = relId[2];
243 if(absId > fEmcCrystals) { // this is CPV
244 fADCchanelCpv = fCalibData->GetADCchannelCpv (module,column,row);
245 fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
246 return fADCpedestalCpv + amp*fADCchanelCpv ;
250 if(absId > fEmcCrystals) // this is CPV
251 return fADCpedestalCpv+ amp*fADCchanelCpv ;
256 //____________________________________________________________________________
257 void AliPHOSClusterizerv1::Exec(Option_t *option)
259 // Steering method to perform clusterization for events
260 // in the range from fFirstEvent to fLastEvent.
261 // This range is optionally set by SetEventRange().
262 // if fLastEvent=-1 (by default), then process events until the end.
264 if(strstr(option,"tim"))
265 gBenchmark->Start("PHOSClusterizer");
267 if(strstr(option,"print")) {
272 GetCalibrationParameters() ;
274 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
276 gime->SetRawDigits(kFALSE);
278 gime->SetRawDigits(kTRUE);
280 if (fLastEvent == -1)
281 fLastEvent = gime->MaxEvent() - 1 ;
283 fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
284 Int_t nEvents = fLastEvent - fFirstEvent + 1;
287 for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
289 gime->Event(ievent ,"D"); // Read digits from simulated data
291 AliRunLoader * rl = AliRunLoader::GetRunLoader(gime->PhosLoader()->GetTitle());
292 rl->GetEvent(ievent);
293 gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
295 fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
299 AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
300 gime->EmcRecPoints()->GetEntries(),ievent));
301 if(AliLog::GetGlobalDebugLevel()>1)
302 gime->EmcRecPoints()->Print();
309 if(strstr(option,"deb"))
310 PrintRecPoints(option) ;
312 //increment the total number of recpoints per run
313 fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
314 fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
317 if(fWrite) //do not unload in "on flight" mode
320 if(strstr(option,"tim")){
321 gBenchmark->Stop("PHOSClusterizer");
322 AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
323 gBenchmark->GetCpuTime("PHOSClusterizer"),
324 gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
328 //____________________________________________________________________________
329 Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
330 Int_t nPar, Float_t * fitparameters) const
332 // Calls TMinuit to fit the energy distribution of a cluster with several maxima
333 // The initial values for fitting procedure are set equal to the positions of local maxima.
334 // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
337 AliPHOSGetter * gime = AliPHOSGetter::Instance();
338 TClonesArray * digits = gime->Digits();
341 gMinuit->mncler(); // Reset Minuit's list of paramters
342 gMinuit->SetPrintLevel(-1) ; // No Printout
343 gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
344 // To set the address of the minimization function
346 TList * toMinuit = new TList();
347 toMinuit->AddAt(emcRP,0) ;
348 toMinuit->AddAt(digits,1) ;
350 gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
352 // filling initial values for fit parameters
353 AliPHOSDigit * digit ;
357 Int_t nDigits = (Int_t) nPar / 3 ;
361 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
363 for(iDigit = 0; iDigit < nDigits; iDigit++){
364 digit = maxAt[iDigit];
369 geom->AbsToRelNumbering(digit->GetId(), relid) ;
370 geom->RelPosInModule(relid, x, z) ;
372 Float_t energy = maxAtEnergy[iDigit] ;
374 gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
377 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
380 gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
383 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
386 gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
389 Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
394 Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
399 gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
400 gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
401 gMinuit->SetMaxIterations(5);
402 gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
404 gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
406 if(ierflg == 4){ // Minimum not found
407 Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
410 for(index = 0; index < nPar; index++){
413 gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index
414 fitparameters[index] = val ;
422 //____________________________________________________________________________
423 void AliPHOSClusterizerv1::GetCalibrationParameters()
425 // Set calibration parameters:
426 // if calibration database exists, they are read from database,
427 // otherwise, they are taken from digitizer.
429 // It is a user responsilibity to open CDB before reconstruction, for example:
430 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
432 AliPHOSGetter * gime = AliPHOSGetter::Instance();
433 // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
435 fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
439 if ( !gime->Digitizer() )
440 gime->LoadDigitizer();
441 AliPHOSDigitizer * dig = gime->Digitizer();
442 fADCchanelEmc = dig->GetEMCchannel() ;
443 fADCpedestalEmc = dig->GetEMCpedestal();
445 fADCchanelCpv = dig->GetCPVchannel() ;
446 fADCpedestalCpv = dig->GetCPVpedestal() ;
450 //____________________________________________________________________________
451 void AliPHOSClusterizerv1::Init()
453 // Make all memory allocations which can not be done in default constructor.
454 // Attach the Clusterizer task to the list of PHOS tasks
456 AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
458 gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
460 AliPHOSGeometry * geom = gime->PHOSGeometry();
462 fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
465 gMinuit = new TMinuit(100);
467 if ( !gime->Clusterizer() ) {
468 gime->PostClusterizer(this);
472 //____________________________________________________________________________
473 void AliPHOSClusterizerv1::InitParameters()
476 fNumberOfCpvClusters = 0 ;
477 fNumberOfEmcClusters = 0 ;
479 fCpvClusteringThreshold = 0.0;
480 fEmcClusteringThreshold = 0.2;
482 fEmcLocMaxCut = 0.03 ;
483 fCpvLocMaxCut = 0.03 ;
491 fEmcTimeGate = 1.e-8 ;
495 fRecPointsInRun = 0 ;
501 SetEventRange(0,-1) ;
503 fIsOldRCUFormat = kFALSE;
506 //____________________________________________________________________________
507 Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
509 // Gives the neighbourness of two digits = 0 are not neighbour but continue searching
511 // = 2 are not neighbour but do not continue searching
512 // neighbours are defined as digits having at least a common vertex
513 // The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
514 // which is compared to a digit (d2) not yet in a cluster
516 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
521 geom->AbsToRelNumbering(d1->GetId(), relid1) ;
524 geom->AbsToRelNumbering(d2->GetId(), relid2) ;
526 if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
527 Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
528 Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
530 if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
531 if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
535 if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
536 rv = 2; // Difference in row numbers is too large to look further
542 if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
549 //____________________________________________________________________________
550 void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
552 // Remove digits with amplitudes below threshold
554 for(Int_t i=0; i<digits->GetEntriesFast(); i++){
555 AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
556 if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
557 (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) )
558 digits->RemoveAt(i) ;
561 for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
562 AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
563 digit->SetIndexInList(i) ;
566 //Overwrite digits tree
567 AliPHOSGetter* gime = AliPHOSGetter::Instance();
568 TTree * treeD = gime->TreeD();
569 treeD->Branch("PHOS", &digits);
571 gime->WriteDigits("OVERWRITE");
572 gime->PhosLoader()->UnloadDigits() ;
574 //____________________________________________________________________________
575 Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
577 // Tells if (true) or not (false) the digit is in a PHOS-EMC module
580 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
582 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
584 if(digit->GetId() <= nEMC ) rv = kTRUE;
589 //____________________________________________________________________________
590 Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
592 // Tells if (true) or not (false) the digit is in a PHOS-CPV module
596 AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
598 Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
600 if(digit->GetId() > nEMC ) rv = kTRUE;
605 //____________________________________________________________________________
606 void AliPHOSClusterizerv1::Unload()
608 AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
609 gime->PhosLoader()->UnloadDigits() ;
610 gime->PhosLoader()->UnloadRecPoints() ;
613 //____________________________________________________________________________
614 void AliPHOSClusterizerv1::WriteRecPoints()
617 // Creates new branches with given title
618 // fills and writes into TreeR.
620 AliPHOSGetter * gime = AliPHOSGetter::Instance();
622 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
623 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
624 TClonesArray * digits = gime->Digits() ;
628 //Evaluate position, dispersion and other RecPoint properties..
629 Int_t nEmc = emcRecPoints->GetEntriesFast();
630 for(index = 0; index < nEmc; index++){
631 AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
632 rp->Purify(fEmcMinE) ;
633 if(rp->GetMultiplicity()>0.) //If this RP is not empty
634 rp->EvalAll(fW0,digits) ;
636 emcRecPoints->RemoveAt(index) ;
640 emcRecPoints->Compress() ;
641 emcRecPoints->Sort() ;
642 // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
643 for(index = 0; index < emcRecPoints->GetEntries(); index++){
644 dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
647 //Now the same for CPV
648 for(index = 0; index < cpvRecPoints->GetEntries(); index++){
649 AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
650 rp->EvalAll(fW0CPV,digits) ;
652 cpvRecPoints->Sort() ;
654 for(index = 0; index < cpvRecPoints->GetEntries(); index++)
655 dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
657 cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
659 if(fWrite){ //We write TreeR
660 TTree * treeR = gime->TreeR();
662 Int_t bufferSize = 32000 ;
663 Int_t splitlevel = 0 ;
666 TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
667 emcBranch->SetTitle(BranchName());
670 TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
671 cpvBranch->SetTitle(BranchName());
676 gime->WriteRecPoints("OVERWRITE");
677 gime->WriteClusterizer("OVERWRITE");
681 //____________________________________________________________________________
682 void AliPHOSClusterizerv1::MakeClusters()
684 // Steering method to construct the clusters stored in a list of Reconstructed Points
685 // A cluster is defined as a list of neighbour digits
688 AliPHOSGetter * gime = AliPHOSGetter::Instance();
690 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
691 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
693 emcRecPoints->Delete() ;
694 cpvRecPoints->Delete() ;
696 TClonesArray * digits = gime->Digits() ;
698 //Remove digits below threshold
699 CleanDigits(digits) ;
701 TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
703 // Clusterization starts
705 TIter nextdigit(digitsC) ;
706 AliPHOSDigit * digit ;
707 Bool_t notremoved = kTRUE ;
709 while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
712 AliPHOSRecPoint * clu = 0 ;
714 TArrayI clusterdigitslist(1500) ;
717 if (( IsInEmc (digit) &&
718 CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
720 CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
721 Int_t iDigitInCluster = 0 ;
723 if ( IsInEmc(digit) ) {
724 // start a new EMC RecPoint
725 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
726 emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
728 emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
729 clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
730 fNumberOfEmcClusters++ ;
731 clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
732 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
734 digitsC->Remove(digit) ;
738 // start a new CPV cluster
739 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
740 cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
742 cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
744 clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
745 fNumberOfCpvClusters++ ;
746 clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
747 clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
749 digitsC->Remove(digit) ;
752 // Here we remove remaining EMC digits, which cannot make a cluster
755 while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
757 digitsC->Remove(digit) ;
761 notremoved = kFALSE ;
768 AliPHOSDigit * digitN ;
770 while (index < iDigitInCluster){ // scan over digits already in cluster
771 digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ;
773 while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
774 Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
776 case 0 : // not a neighbour
778 case 1 : // are neighbours
779 if (IsInEmc (digitN))
780 clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
782 clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
784 clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
786 digitsC->Remove(digitN) ;
788 case 2 : // too far from each other
797 } // loop over cluster
808 //____________________________________________________________________________
809 void AliPHOSClusterizerv1::MakeUnfolding()
811 // Unfolds clusters using the shape of an ElectroMagnetic shower
812 // Performs unfolding of all EMC/CPV clusters
814 AliPHOSGetter * gime = AliPHOSGetter::Instance();
816 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
818 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
819 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
820 TClonesArray * digits = gime->Digits() ;
822 // Unfold first EMC clusters
823 if(fNumberOfEmcClusters > 0){
825 Int_t nModulesToUnfold = geom->GetNModules() ;
827 Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
829 for(index = 0 ; index < numberofNotUnfolded ; index++){
831 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
832 if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
835 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
836 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
837 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
838 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
840 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
841 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
842 emcRecPoints->Remove(emcRecPoint);
843 emcRecPoints->Compress() ;
845 fNumberOfEmcClusters -- ;
846 numberofNotUnfolded-- ;
849 emcRecPoint->SetNExMax(1) ; //Only one local maximum
853 delete[] maxAtEnergy ;
856 // Unfolding of EMC clusters finished
859 // Unfold now CPV clusters
860 if(fNumberOfCpvClusters > 0){
862 Int_t nModulesToUnfold = geom->GetNModules() ;
864 Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
866 for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
868 AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
870 if(recPoint->GetPHOSMod()> nModulesToUnfold)
873 AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint*>(recPoint) ;
875 Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
876 AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
877 Float_t * maxAtEnergy = new Float_t[nMultipl] ;
878 Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
880 if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
881 UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
882 cpvRecPoints->Remove(emcRecPoint);
883 cpvRecPoints->Compress() ;
885 numberofCpvNotUnfolded-- ;
886 fNumberOfCpvClusters-- ;
890 delete[] maxAtEnergy ;
893 //Unfolding of Cpv clusters finished
897 //____________________________________________________________________________
898 Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r)
900 // Shape of the shower (see PHOS TDR)
901 // If you change this function, change also the gradient evaluation in ChiSquare()
903 Double_t r4 = r*r*r*r ;
904 Double_t r295 = TMath::Power(r, 2.95) ;
905 Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
909 //____________________________________________________________________________
910 void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
912 AliPHOSDigit ** maxAt,
913 Float_t * maxAtEnergy)
915 // Performs the unfolding of a cluster with nMax overlapping showers
917 AliPHOSGetter * gime = AliPHOSGetter::Instance();
919 const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
921 const TClonesArray * digits = gime->Digits() ;
922 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
923 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
925 Int_t nPar = 3 * nMax ;
926 Float_t * fitparameters = new Float_t[nPar] ;
928 Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
930 // Fit failed, return and remove cluster
931 iniEmc->SetNExMax(-1) ;
932 delete[] fitparameters ;
936 // create ufolded rec points and fill them with new energy lists
937 // First calculate energy deposited in each sell in accordance with fit (without fluctuations): efit[]
938 // and later correct this number in acordance with actual energy deposition
940 Int_t nDigits = iniEmc->GetMultiplicity() ;
941 Float_t * efit = new Float_t[nDigits] ;
942 Float_t xDigit=0.,zDigit=0.,distance=0. ;
943 Float_t xpar=0.,zpar=0.,epar=0. ;
945 AliPHOSDigit * digit = 0 ;
946 Int_t * emcDigits = iniEmc->GetDigitsList() ;
950 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
951 digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
952 geom->AbsToRelNumbering(digit->GetId(), relid) ;
953 geom->RelPosInModule(relid, xDigit, zDigit) ;
957 while(iparam < nPar ){
958 xpar = fitparameters[iparam] ;
959 zpar = fitparameters[iparam+1] ;
960 epar = fitparameters[iparam+2] ;
962 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
963 distance = TMath::Sqrt(distance) ;
964 efit[iDigit] += epar * ShowerShape(distance) ;
969 // Now create new RecPoints and fill energy lists with efit corrected to fluctuations
970 // so that energy deposited in each cell is distributed betwin new clusters proportionally
971 // to its contribution to efit
973 Float_t * emcEnergies = iniEmc->GetEnergiesList() ;
977 while(iparam < nPar ){
978 xpar = fitparameters[iparam] ;
979 zpar = fitparameters[iparam+1] ;
980 epar = fitparameters[iparam+2] ;
983 AliPHOSEmcRecPoint * emcRP = 0 ;
985 if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
987 if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
988 emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
990 (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
991 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
992 fNumberOfEmcClusters++ ;
993 emcRP->SetNExMax((Int_t)nPar/3) ;
995 else{//create new entries in fCpvRecPoints
996 if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
997 cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
999 (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
1000 emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
1001 fNumberOfCpvClusters++ ;
1005 for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
1006 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
1007 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1008 geom->RelPosInModule(relid, xDigit, zDigit) ;
1009 distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
1010 distance = TMath::Sqrt(distance) ;
1011 ratio = epar * ShowerShape(distance) / efit[iDigit] ;
1012 eDigit = emcEnergies[iDigit] * ratio ;
1013 emcRP->AddDigit( *digit, eDigit ) ;
1017 delete[] fitparameters ;
1022 //_____________________________________________________________________________
1023 void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag)
1025 // Calculates the Chi square for the cluster unfolding minimization
1026 // Number of parameters, Gradient, Chi squared, parameters, what to do
1028 TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
1030 AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
1031 TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
1035 // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
1037 Int_t * emcDigits = emcRP->GetDigitsList() ;
1039 Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
1041 Float_t * emcEnergies = emcRP->GetEnergiesList() ;
1043 const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
1048 for(iparam = 0 ; iparam < nPar ; iparam++)
1049 Grad[iparam] = 0 ; // Will evaluate gradient
1053 AliPHOSDigit * digit ;
1056 for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
1058 digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
1064 geom->AbsToRelNumbering(digit->GetId(), relid) ;
1066 geom->RelPosInModule(relid, xDigit, zDigit) ;
1068 if(iflag == 2){ // calculate gradient
1071 while(iParam < nPar ){
1072 Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
1074 distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
1075 distance = TMath::Sqrt( distance ) ;
1077 efit += x[iParam] * ShowerShape(distance) ;
1080 Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
1082 while(iParam < nPar ){
1083 Double_t xpar = x[iParam] ;
1084 Double_t zpar = x[iParam+1] ;
1085 Double_t epar = x[iParam+2] ;
1086 Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
1087 Double_t shape = sum * ShowerShape(dr) ;
1088 Double_t r4 = dr*dr*dr*dr ;
1089 Double_t r295 = TMath::Power(dr,2.95) ;
1090 Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
1091 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
1093 Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
1095 Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
1097 Grad[iParam] += shape ; // Derivative over energy
1104 while(iparam < nPar ){
1105 Double_t xpar = x[iparam] ;
1106 Double_t zpar = x[iparam+1] ;
1107 Double_t epar = x[iparam+2] ;
1109 Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
1110 distance = TMath::Sqrt(distance) ;
1111 efit += epar * ShowerShape(distance) ;
1114 fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
1115 // Here we assume, that sigma = sqrt(E)
1120 //____________________________________________________________________________
1121 void AliPHOSClusterizerv1::Print(const Option_t *)const
1123 // Print clusterizer parameters
1126 TString taskName(GetName()) ;
1127 taskName.ReplaceAll(Version(), "") ;
1129 if( strcmp(GetName(), "") !=0 ) {
1131 message = "\n--------------- %s %s -----------\n" ;
1132 message += "Clusterizing digits from the file: %s\n" ;
1133 message += " Branch: %s\n" ;
1134 message += " EMC Clustering threshold = %f\n" ;
1135 message += " EMC Local Maximum cut = %f\n" ;
1136 message += " EMC Logarothmic weight = %f\n" ;
1137 message += " CPV Clustering threshold = %f\n" ;
1138 message += " CPV Local Maximum cut = %f\n" ;
1139 message += " CPV Logarothmic weight = %f\n" ;
1141 message += " Unfolding on\n" ;
1143 message += " Unfolding off\n" ;
1145 message += "------------------------------------------------------------------" ;
1148 message = " AliPHOSClusterizerv1 not initialized " ;
1150 AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
1155 fEmcClusteringThreshold,
1158 fCpvClusteringThreshold,
1164 //____________________________________________________________________________
1165 void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
1167 // Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
1169 AliPHOSGetter * gime = AliPHOSGetter::Instance();
1171 TObjArray * emcRecPoints = gime->EmcRecPoints() ;
1172 TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
1174 AliInfo(Form("\nevent %d \n Found %d EMC RecPoints and %d CPV RecPoints",
1175 gAlice->GetEvNumber(),
1176 emcRecPoints->GetEntriesFast(),
1177 cpvRecPoints->GetEntriesFast() )) ;
1179 fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
1180 fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
1183 if(strstr(option,"all")) {
1184 printf("\n EMC clusters \n") ;
1185 printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
1187 for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
1188 AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
1190 rp->GetLocalPosition(locpos);
1192 rp->GetElipsAxis(lambda);
1195 primaries = rp->GetPrimaries(nprimaries);
1196 printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
1197 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
1198 locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
1200 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
1201 printf("%d ", primaries[iprimary] ) ;
1206 //Now plot CPV recPoints
1207 printf("\n CPV clusters \n") ;
1208 printf("Index Ene(MeV) Module X Y Z \n") ;
1209 for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
1210 AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
1213 rp->GetLocalPosition(locpos);
1215 printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
1216 rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
1217 locpos.X(), locpos.Y(), locpos.Z()) ;